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Abstract 



This work reports on the discovery of HESS J 15 14— 591, a VHE y-ray source found at the pulsar 
wind nebula (PWN) MSH 15-52 and its associated pulsar PSR B1509-58. The discovery was 
made with the High Energy Stereoscopic System (H.E.S.S.), which currently provides the most 
sensitive measurement in the energy range of about 0.2-100 TeV. This analysis is the first to 
include all H.E.S.S. data from observations dedicated to MSH 15—52. The data was taken 
in 2004 from March 26 to July 20, with a total live-time of 26.14h. The y-ray signal was 
detected with a statistical significance of 32 standard deviations. The intensity distribution 
shows an elliptical extension with the major axis oriented in a southeast direction. The standard 
deviations of a Gaussian fit function are 6.5' ± 0.55(j^( ± 0. l[y^^ and 2.3' ± 0A[^^^ ± 0. 1 ^y^j for the 
major- and minor axis, respectively. The y-ray emission extends in direction of the pulsar jet, 
previously resolved in X-rays. This becomes more apparent after image deconvolution. The 
emission region along the jet axis decreases with increasing energy. The corresponding flux 
above ITeV is (4.4 ± 0.2stat ± 1-Osyst) x 10~^^cm~^s~^ The energy spectrum obeys a power 
law with a differential flux at 1 TeV of (5.8 ± 0.2stat ± 1.3syst) x 10-12^^-25-1 jgy^i ^nd a 
photon index of 2.32 ± 0.04stat ± 0. lOsyst- The y-ray light curve with periodicity according to 
PSR B1509— 58 yields a uniform distribution. An upper limit of 11. Ox 10 ^^cm^^s ^ for the 
pulsed y-ray flux from PSR B 1509— 58 was calculated with a confidence level of 99%. 

In addition to these results the following subjects are discussed: previous observations of 
MSH 15 52 and PSR B 1509 58, the theory of pulsars, PWNs and their y-ray production, the 
imaging air Cherenkov technique for the detection of y radiation in the earth's atmosphere, the 
H.E.S.S. experiment and its data analysis, the first (Richardson-Lucy) deconvolution of VHE 
gamma-ray maps, the analysis of H.E.S.S. data for pulsed emission from pulsars using radio 
ephemeris. 

The results are discussed within the framework of PWNs and are explained by inverse 
Compton scattering of leptons. A hadronic component in MSH 15—52 is not excluded, but 
its y-ray emission would not be significant. Moreover, it is concluded that advection is the 
dominant transport mechanism over diffusion in the magnetized flow of the pulsar wind from 
PSR B 1509— 58. A correlation analysis with the Chandra X-ray data suggests that the y radia- 
tion is emitted from the region of PSR B 1509— 58, but not from the neighboring optical nebula 
RCW 89. 
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Chapter 1 
Introduction 



"We owe our existence to stars, because they make the atoms of which we are 
formed. So if you are romantic you can say we are literally starstuff. If you 're less 
romantic you can say we're the nuclear waste from the fuel that makes stars shine. 

We 've made so many advances in our understanding. A few centuries ago, the 
pioneer navigators learnt the size and shape of our Earth, and the layout of the 
continents. We are now just learning the dimensions and ingredients of our entire 
cosmos, and can at last make some sense of our cosmic habitat." 

— Sir Martin Rees, British astrophysicist and president of the Royal Society 



Astronomy is and always has been a central discipline of natural science, driven by fundamen- 
tal questions and exiting answers. The first recorded astronomical achievements date back to 
early cultures such as the Babylonians, Egyptians and Chinese. Further progress was made in 
the Renaissance, when the heliocentric model of the solar system was proposed by Nicolaus 
Copernicus, Galileo Galilei and Johannes Kepler. The use of the telescope for astronomic ob- 
servations by Galilei marks the beginning of experimental astronomy. Since then, astronomy 
has evolved rapidly. For example, the introduction of spectroscopy and photography by Joseph 
Fraunhofer in 1814 laid the foundations for a "New Astronomy" and astrophysics by provid- 
ing the means for determining the chemical composition of astronomical objects. Moreover, it 
paved the way for the determination of red shifts by Vesto Slipher in 1912, which allowed for 
such far-reaching conclusions as the expansion of the universe by Hubble in 1929. An astro- 
physical revolution began in the second half of the 20th century, when new types of telescopes 
became available, owed to technological advances, with which the full range of the electro- 
magnetic spectrum could be explored. Radio telescopes permitted the discovery of the cosmic 
microwave background radiation in 1965 and pulsars in 1967, both of which were honored with 
the Nobel prize; infrared telescopes revealed the view through vast dust clouds to previously 
hidden objects; X- and 7-ray satellites provided pictures of the non-thermal universe and its 
most violent processes, such as 7-ray bursts first observed in 1967, active galactic nuclei or 
super nova remnants; also new experiments in the rising field of astroparticle physics provided 
fresh insights into the non-thermal universe, e.g. Kamiocande, the Irvine-Michigan-Brookhaven 
detector and the scintillator experiment at Baksan by the detection of neutrinos from the super- 
nova SN 1987 A. Although these examples only name some of astrophysical milestones, they 
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demonstrate that the exploration of new fields of astronomy can lead to outstanding discoveries 
with great impact on the understanding of the universe. 

In this respect, imaging atmospheric Cherenkov telescopes (lACT) have been developed 
for the exploration of the very high energy (VHE) y-ray sky which extends from lOGeV to 
lOOTeV. Therefore, lACTs currently provide a window to the highest available y-ray energies. 
Since this / radiation is produced where highly accelerated particles interact with their environ- 
ment, TeV 7 radiation provides important information about the acceleration mechanisms for 
the primary particles. In comparison, 7-ray astronomy using lACTs is a relatively young field 
which achieved its breakthrough in 1989 with the discovery of TeV 7 radiation from the Crab 
Nebula by the Whipple collaboration. Since then, lACT based 7 astronomy has progressed sig- 
nificantly. While it took several weeks for the first detection of the Crab Nebula which is known 
as the strongest TeV 7-ray source, nowadays lACTs can detect the Crab Nebula within less than 
a minute. Moreover, while previously only a handful of TeV sources could be detected, today 
about 50 TeV 7-ray sources have been established and almost monthly the detection of a new 
source is reported. Many of these new detections are owed to the High Energy Stereoscopic 
System (H.E.S.S.), which is currently one of the most sensitive LACTs worldwide. During its 
first few years of operation, starting in 2002, H.E.S.S. has detected or confirmed more than 40 
sources of TeV 7 radiation (cf. Hofmann [2005] and Aharonian et al. [2005c]). 

One of these sources is the pulsar wind nebula (PWN) MSH 15—52, at a distance of about 
5.2 kpc from earth. PWNs, also called Plerions, are very unique but also very rare objects of 
which only about 50 have been identified, all within the Galaxy. The central object in a PWN 
is a pulsar which exposes extreme condition to its environment and generates 7 radiation in its 
vicinity, in particular by the emission of a wind of VHE particles. The pulsar associated with 
MSH 15—52 is PSR B 1509— 58, which is one of the most energetic pulsars known. Many ob- 
servations of MSH 15—52 have been conducted from radio to 7-ray energies since its discovery 
in 1961. They have shed light on the violent emission processes. However, the observations at 
the highest energies in the TeV range, which are crucial for the understanding of the processes 
in MSH 15—52 and in PWNs in general, were rather incomplete, since no experiment had been 
able to provide them. Therefore, when H.E.S.S. came into operation with its unprecedented 
sensitivity, MSH 15—52 was scheduled as one of the first targets for thorough observation over 
a period of several month. This work presents this data from the first H.E.S.S. studies. It dis- 
cusses details of the analysis, the detection, and further results as well as possible implications 
for the astrophysical processes involved. As it turns out, MSH 15—52 is one of the strongest 
TeV 7-ray sources that have ever been detected. For a comprehensive discussion, also an intro- 
duction to the system of MSH 15—52 and PSR B1509— 58, the imaging atmospheric Cherenkov 
technique and the H.E.S.S. experiment is given in advance in the following chapters. 
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Chapter 2 

MSH 15-52 and PSR B1509-58 



2.1 Review of Previous Observations 

The supernova remnant (SNR) MSH 15—52, also known as G320.4— 1.2, was first discovered 
by Mills et al. [1961] as an extended radio source. Later, radio observations by Caswell et al. 
[1981] did resolve individual, non-thermal radio features of the SNR, one of them coinciding 
with the optical nebula RCW 89 in the northwest. X-ray observations in the energy range 
from 0.2-4 keV with the Einstein satellite by Seward and Harnden [1982] revealed a possibly 
associated pulsar with an increasing period of about 150 ms, a characteristic age of about 1 .6 kyr 
and a surrounding nebula. Calculations showed that the nebula could easily be powered by the 
pulsar. Subsequent radio observations by Manchester et al. [1982] confirmed the existence of 
the pulsar PSR B 1509-58. Nevertheless, it was still unclear whether MSH 15-52 and PSR 
B 1509— 58 are associated and whether MSH 15—52 was a pulsar wind nebula (PWN) similar 
to the Crab Nebula. A main argument against an association was the unexplained difference 
between the apparent age of the pulsar and the supernova remnant. However, these questions 
have been resolved in later years in favor of an association. In 1983 Seward et al. [1983] reported 
on infrared spectral lines in the region of RCW 89 and provided an optical image in which RCW 
89 and a new northerly filament became apparent. The detection of TeV y radiation from the 
region of MSH 15—52 was reported by Sako et al. [2000]. Caraveo et al. [1994] also suggested 
an optical counterpart of PSR B 1509— 58. However, Kaplan and Moon [2006] found a more 
likely infrared counterpart which was hidden by the object proposed by Caraveo et al. [1994]. 

A historical association of MSH 1 5 —52 with the ^uest star of AD 1 85 December 7, which 
was witnessed by Chinese astronomers according to the Houhanshu, has been pointed out by 
Strom [1994] and Schaefer [1995]. The age and position of the guest star would allow for an 
interpretation of it as the supernova of MSH 15—52. 

MSH 15—52 is located in the galactic plane with an offset of 40° from the galactic center. 
Fig. 2.1 shows this location overlaid to the galactic map of the Third EGRET Catalog, which 
was obtained by the Energetic Gamma Ray Experiment Telescope (EGRET) and which covers 
the energy range from 100 MeV to 30GeV. MSH 15—52 does not coincide with any of the 
sources shown. The estimated distance of MSH 15—52 from earth is 5.2 kpc. 

Tbl. 2.1 and 2.2 summarize basic parameters of MSH 15-52 and PSR B1509-58 which 
have been determined from observations at different wavelengths discussed below. 

^^The ancient Chinese term for a star that newly appears and is visible for a short time. 
^^The Houhanshu are Chinese records of the Later Han dynasty. 
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Figure 2.1: Galactic map of the Third EGRET Catalog overlaid with the location of MSH 15—52. 
(Figure taken from NASA [2006a].) 

Table 2.1: Parameters of the supernova remnant MSH 15— 52 as found by Gaensler et al. [1999]. 



Parameter 


Variable 


Value 


Celestial coordinates 


RA (J2000) 






Dec (J2000) 


-59°16f3 


Galactic coordinates 


/ 


320.31° 




b 


-1.31° 


Diameter 


D 


5' 


Age 


Tn 


(2-20) xlO^ yr 


Distance 


d 


(5.2±1.4)kpc 


IVIagnetic field of the PWN 


BpWN 


-5-8 AiG 



Table 2.2: Parameters of PSR B 1509-58 for epoch 48355.0000 MJD determined by Kaspi et al. 
[1994] from radio observations at 843 MHz with the Molonglo Observatory Synthesis Telescope 
(MOST, Mills [1981]). A * marks parameters found by Gaensler et al. [2002]. 



Parameter 


Variable 


Value 


Celestial coordinates 


RA (J2000) 


15'^13'"55.^62±0.W 




Dec (J2000) 


-59°08'09."0±lf'0 


Period 


P 


151 ms 


Frequency 


f 


6.6375697328(8) s-i 


First Frequency derivative 


f 


-6.7695374(4) x lO-^^s 


Second Frequency derivative 


f 


1.9587(9) X 10-2is-3 


Third Frequency derivative 


f 


-1.02(25) X 10-3is-4 


Breaking index 


n 


2.837 ±0.001 


Second deceleration parameter 


m 


14.5 ±3.6 


Dispersion measure 




(253.2 ± 1.9) pc cm-3 


Characteristic age 


Tc 


1700 yr 


Spin-down luminosity* 


E 


1.8 X 10^7 erg s^i 


Surface magnetic field* 


Bp 


1.5x1013g 
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2.1.1 MSH 15-52 
Radio Observations 

The radio observations of MSH 15—52 show an unusual appearance of two separated regions 
with non-thermal emission. The situation is illustrated in Fig. 2.2, which shows the radio map 
from Whiteoak and Green [1996] overlaid with X-ray contours from Trussoni et al. [1996]. 
While the southern region approximates a partial shell, the northern region is rather compact 
and coincides with the optical nebula RCW 89, which contains an unusual ring of radio clumps. 
From the region of PSR B 1509— 58 no signal is seen. With standard parameters for the super- 
nova and the interstellar medium, the size of the SNR suggests an age of ~6-20kyr, which is 
an order of magnitude larger than the spin-down age of the pulsar of 1 .7 kyr and which is in 
contradiction to the association of MSH 15—52 and PSR B 1509— 58. However, Gaensler et al. 
[1999] have confirmed by Hi absorption measurements that the observed components are 
part of a single SNR. 



Figure 2.2: Radio and X-ray image of MSH 
15—52. The gray scale corresponds to the 
MOST radio observation at 843 MHz by 
Whiteoak and Green [1996]. The white con- 
tour lines at the levels of 0.5, 1, 1.5, 2, 5, 
10, 20, 30 and 40 in arbitrary units represent 
the ROSAT observations from Trussoni et al. 
[1996]. The position of PSR B 1509-58 is 
marked by the cross. The black box marks 
the Chandra ACTS -I field of view. (Figure 
taken from Gaensler et al. [2002].) 



15^16™ 15™ 14™ 13™ 12' 

RA (J2000) 



Infrared Observations 

Seward et al. [1983] have detected infrared spectral lines of many elements from the region of 
RCW 89. Among them Fell at 1.76 /im, which was detected for the first time from an SNR. 
The spectrum is clearly non-thermal and typical for a reddened, high density and coUisionally 
excited nebula at a distance of 5 kpc. 

Optical Observations 

Fig. 2.3 from Seward et al. [1983] shows the area of MSH 15—52 from a 90min R-band 
(A =610-700 nm) plate overlaid with the X-ray contours from the Einstein satellite. The re- 
gion of RCW 89 and a filament 9' to the northeast are apparent. The filament coincides with 

^^H I and H n denote neutral and ionized hydrogen, respectively. H i measurements look for the 21 cm hne of 
the forbidden hyper fine transition. It is used to determine the density and velocity of hydrogen. 
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the radio arc of MSH 15—52, which extends eastwards from RCW 89. Due to the good agree- 
ment with the radio observations, Seward et al. [1983] concluded that the whole northwest arc 
is a single entity. They also pointed out that the association of the SNR and PSR B 1509— 58 
can still be explained if a local low density of the interstellar medium is assumed, in which the 
SNR could have expanded unusually rapidly. The low density could be explained by an earlier 
supernova in the same region. 




Figure 2.3: Optical image of the region of MSH 
15-52 from a 90 min R-band (A =610-700 nm) 
plate overlaid with X-ray contours from the Ein- 
stein satellite. The region of RCW 89 and a fil- 
ament 9' to the northeast are apparent. (Image 
taken from Seward et al. [1983].) 



X-ray Observations 

The X-ray emission from MSH 15—52 has several components. The component in the region of 
PSR B 1509— 58 shows highest intensity and pulsed emission. Another component is elongated 
and extending southeast of the pulsar. They have been identified as a PWN and as a jet by 
several authors including Gaensler et al. [2002] and Forot et al. [2006]. A third component 
is in the region of RCW 89, which also reaches a considerable peak intensity, but in contrast 
to the former has a thermal spectrum. There is also a faint diffuse X-ray emission throughout 
the whole region of the SNR, which can be considered as a fourth component. According 
to Trussoni et al. [1996], the spectrum is compatible with thermal and non-thermal emission 
and could therefore correspond to thermal emission from the SNR blast wave (Gaensler et al. 
[2002]). 

The four Chandra satellite images of different energy bands in Fig. 2.4 from Gaensler et al. 
[2002] provide detailed insights into the structure of the PWN and the energy spectrum of the 
individual components. Clearly, the pulsar has the hardest spectrum, followed by the jet and 
the PWN. The region of RCW 89 and the contained X-ray and radio clumps disappear with 
increasing energy. A closer analysis by Gaensler et al. [2002] showed that the spectrum in 
this region is dominated by emission lines, in contrast to the spectrum in the pulsar region. 
The energy spectrum of the diffuse PWN was determined with a photon index of 2.05 ±0.04. 
Fig. 2.5 shows the same data in a smaller region surrounding PSR B 1509— 58 in the energy 
range from 0.3-8.0 keV. It reveals several distinct features which are discussed in Gaensler et al. 
[2002]. Feature A corresponds to the pulsar. Noteworthy is feature 5 which refers to a faint 
circular arc of emission. It is approximately centered at the pulsar, with a radius of 17" and an 
angle subtended at the pulsar of ~ 1 10°. 

Moreover Gaensler et al. [2002] concluded that the PWN has a magnetic field of Bp^^ ~ 
8 iiG and that the outflow of the jet has a velocity of 0.2c, carrying away at least 0.5% of the 
pulsar's spin-down luminosity. 
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Figure 2.4: Chandra X-ray images of MSH 15-52 
showing different energy bands. The images are expo- 
sure corrected, convolved with a Gaussian of FWHM 
of 10" and displayed using a linear transfer function. 
(Figure taken from Gaensler et al. [2002].) 




RA (J2000) 



Figure 2.5: High resolution Chandra X-ray 
image of MSH 15—52 showing the immedi- 
ate surrounding of PSR B 1509-58 (corre- 
sponding to feature A) in the energy range 
0.3-8.0 keV. The image is exposure corrected 
but no smoothing was appUed. The labeled 
features are discussed in the text. (Figure 
taken from Gaensler et al. [2002].) 



Further studies of the Chandra and ROSAT data led DeLaney et al. [2006] to the discovery 
that bright knots within 20" of the pulsar show an even higher outflow velocity, up to 0.6c. Also, 
significant time variability of the brightness has been found. For example, the brightness of the 
jet has increased by 30% within 9 years. 

Recent observations with the INTEGRAL X-ray satellite in 2005 in the energy range from 
20-200 keV reported by Forot et al. [2006] predict a similar high outflow velocity of 0.3-0.5c 
corresponding to parent electron energies of 400-730 TeV and a mean magnetic field strength 
of 22-33 juG with a systematic uncertainty of the later of 27%. Comparisons with other X-ray 
measurements show a jet length (L) scaling with the X-ray energy (E) as L<^ E^. The analysis 
also shows the source extension with a standard deviation of the major axis of 5'.53 ± O'.Ol. The 
energy spectrum of the PWN was found to obey a power law with a photon index of 2. 12 ± 0.05 
in agreement with the Chandra observations. However, in contrast to previous measurements 
by the BeppoSAX satellite the INTEGRAL data suggests a spectral break near 160 keV. Also 
pulsed emission from PSR B 1509— 58 was clearly resolved. 



TeV 7-Ray Observations 

/radiation from MSH 15—52 was first predicted by Du Plessis et al. [1995] based on similarities 
to the Crab Nebula. These were the first prediction of VHE 7 radiation from a PWN other 
than the Crab nebula. Sako et al. [2000] reported the first evidence for 7 radiation from MSH 
15—52 based on observations with the lACT CANGAROO. An excess with a significance of 4.1 
standard deviations was found in the region of PSR B 1509— 58 (Fig. 2.6). The corresponding 
flux above 1.9 TeV was determined to (2.9 ±0.7) x 10~^^cm"^s^^ A magnetic field strength 
of 5juG was estimated. No pulsed 7-ray emission was found. 
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Figure 2.6: y-Ray map of MSH 15-52 

showing the first evidence for TeV 7 radi- 
ation from this region. (Figure taken from 
Sako et al. [2000].) 



Right Ascension (degree) 



2.1.2 PSR B1509-58 

PSR B 1509— 58 is one of the most energetic pulsar known, having very high spin-down lumi- 
nosity and magnetic fields. It is also one of the youngest pulsars, and for a young pulsar its spin 
is extraordinarily stable. No glitches have occurred within 24 years since its discovery in 1982. 
Owing to this stability, Kaspi et al. [1994] have determined the spin parameters with very high 
precision from radio observations over six years (Tbl.2.2). The second deceleration parameter 
m is consistent with a constant breaking index and magnetic moment on timescales of kyr. 

PSR B1509— 58 has been detected from radio to y-ray energies i.e. in the range from 10~^^ 
to 10^ eV. Fig. 2.7 from Thompson et al. [1999] shows the spectral energy distribution in com- 
parison to other known 7-ray pulsars. Characteristic for PSR B 1509— 58 is the low optical and 
7-ray flux indicated by the upper limits in Fig. 2.7. The light curve of PSR B 1509— 58 is shown 
in Fig. 2.8 from radio to low 7-ray energies in comparison to light curves of the Crab and Vela 
Pulsar. A phase lag is apparent which increases with energy. A detailed investigation of the 
light curve of PSR B 1509— 58 near its cutoff energy was done by Kuiper et al. [1999]. With an 
analysis of CGRO data from COMPTEL and EGRET, they detected a signal with a modu- 
lation significance of 5.6 standard deviations in the energy range from 0.75-30 MeV and found 
the cutoff energy at ~ 10 MeV. They also found indications for a double-peaked profile at X-ray 
energies. Fig. 2.9 shows the light curve of this analysis in the energy range from 0.75-10 MeV 
along with the X-ray light curves. The phase lag of about 0.3 to the radio phase is also indicated. 

According to Harding et al. [1997], an explanation for the low 7 radiation at high energies 
could be the extraordinary strong magnetic field, which causes a spectral cutoff at a few MeV 
due to photon splitting. 

The infrared counterpart of PSR B 1509—58 is faint but was probably discovered recently 
by Kaplan and Moon [2006]. Fig. 2.10 shows the Kj-band (1.99-2. 30 /im) of these observa- 
tions with the Persson's Auxiliary Nasmyth Infrared Camera (PANIC, Martini et al. [2004]). 
The source labeled with "A" is the proposed counterpart with a magnitude of ~ 19.4, which 

"♦^CGRO is the acronym for the Compton Gamma Ray Observatory, which was a satellite mission from 1991 to 
2000. Two of its instruments were the Imaging Compton Telescope (COMPTEL) and the Energetic Gamma Ray 
Experiment Telescope (EGRET). 
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Figure 2.7: Spectral energy 
distribution of PSR B 1509-58 
from multiwavelength obser- 
vation in comparison to other 
pulsars. Note that for PSR 
B 1509-58 the optical data 
point and the data points at 
energies > 1 MeV only indicate 
upper limits. (Figure taken 
from Thompson et al. [1999].) 
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Figure 2.8: Light curve of PSR 
B1509— 58 in comparison to the 
hght curves of the Crab and the 
Vela Pulsar. Characteristic for 
PSR B 1509-58 is a phase lag 
increasing with energy and a 
low optical emission. (Figure 
taken from NASA [2006b].) 
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Figure 2.9: Light curves of 
PSR B 1509-58 at high ener- 
gies: 0.7-2.2 keV (ASCA, Saito 
et al. [1997]), 2-16 keV (RXTE, 
Rots et al. [1998]), >32keV 
(BATSE, Rots et al. [1998]) 
and 0.75-10 MeV (COMTEL, 
Kuiper et al. [1999]). (Figure 
taken from Kuiper et al. [1999].) 



coincides with the pulsar's X-ray position determined by Gaensler et al. [2002]. Kaplan and 
IVIoon [2006] have also pointed out that the previously proposed optical counterpart by Caraveo 
et al. [1994] is likely to be the object labeled "C]VIB94" and not PSR B 1509-58. The optical 
measurement should therefore rather be considered as an upper limit. So the detection of an 
optical counterpart remains a task for future observations. 

Further pulsar parameters derived from different observations are presented in Tbl. 2.2. 



Figure 2.10: Potential infrared counterpart 
of PSR B 1509-58, visible with a magni- 
tude of ~ 19.4 on a K^-band (1.99-2.30 urn) 
image from the PANIC infrared camera re- 
ported by Kaplan and Moon [2006]. Source 
"A" is the proposed counterpart. The object 
labeled "CMB94" is probably the optical ob- 
ject previously discovered by Caraveo et al. 
[1994]. 
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2.2 Model of the Pulsar 

To understand how the pulsar parameters of Tbl. 2.2 can be derived from observations, it is 
necessary to be familiar with pulsar models. The basic concepts are discussed in this section. 
Further details can be found in a book by Lyne and Graham-Smith [1998]. 

The models assume that pulsars are rotating neutron stars with high magnetic fields. The 
rotation of the magnetic fields and the supply of particles from the neutron stars' surface can 
eventually lead to the periodic emission of strong electromagnetic radiation analogous to a 
lighthouse. 

2.2.1 Formation and Inner Structure 

Although only about a dozen out of more than 700 known pulsars appear to be convincingly 
associated directly with supernova remnants, neutron stars are commonly believed to originate 
in a supernova from the collapse of a star's core. Thereby the mass of the core determines its 
density and whether the core forms a white dwarf, a neutron star or a black hole. An important 
criterion is the Chandrasekhar limit (Mch), which is given by 

fhc\^ 1 

^Ch^[^j ^ = 1.44M0, (2.1) 

where h is the reduced Planck constant, c is the speed of light, G is the gravitational constant, 
is the mass of a proton and Mq = 2 x 10^^ kg is the mass of the sun. For masses exceeding Mch, 
the gravitational pressure exceeds the electron degeneracy pressure^^ inside the core, such that 
the atoms are compressed. Electron capture by the nuclei is the consequence that eventually 
leads to the formation of neutrons by inverse j8 decay (p+ + e~ — > n + Ve), resulting in an 
extremely compact state of matter — a neutron superfluid. With an increasing mass of the core, 
the ratio of neutrons to atoms also increases, while the radius decreases. However, an upper 
mass limit for a neutron star is reached at ~ 2.5Mq, where gravitation compresses the neutron 
star below its Schwarzschild radius, converting it in a black hole. So the radius of a neutron star 
ranges from between 10 to 20 km only. The moment of inertia (/) is therefore of the order of 
3 X 10^^ g/cm^, which approximately equals the moment of inertia of the earth. 

The assumed structure of a neutron star is shown in Fig. 2. 1 1 . It consists of different layers. 
The outer layer, i.e. the crust, presumably contains a rigid lattice of heavy atoms, such as iron 
with a density of ~ 10^^ g/cm?. The middle layer contains a neutron superfluid with a density 
of ~ lO^'* g/cm^. The inner region might contain a solid core with a density up to 10^^ g/cm^. 

2.2.2 Spin-down Luminosity 

The conservation of the angular momentum and magnetic flux of the progenitor star can lead to 

high angular velocities {Q. = 27tf), i.e. short rotational periods {P = 1//) and extremely high 
magnetic fields (5) of the pulsar. Typical values range from a few milliseconds to a few seconds 
for P and from 10^ to 10^^ Gauss for B. Such magnetic fields are extremely high. Their energy 
density is equivalent to a mass density of 1 kg/cm? (Ruderman [1974]). 

^^The electron degeneracy pressure is caused by the Pauli exclusion principle, which states that two electrons 
cannot occupy the same quantum state at the same time. 
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Figure 2.11: Inner structure of a neutron star. 
(Figure taken from NASA/HEASARC [2006].) 



By the rotation of their strong magnetic fields, pulsars lose a significant amount of their 
rotational energy {Erot = 5^^^) by electromagnetic dipole radiation. The total loss of rotational 
energy is called spin-down luminosity (E) and is given by 

d . P 

E = - —Erot = -I^^ = ■ (2.2) 

Thus, it can be determined by a measurement of Q. and Cl. For a magnetic field with the moment 
(M^) perpendicular to the axis of rotation, the electromagnetic radiation is |M^c~^^2^ and the 
relation between Cl and D. reads 

ma = fa] (2.3) 

3 c3 ^ 



2.2.3 Characteristic Age and Breaking Index 

Observations have shown that pulsars dissipate their energy not exclusively by electromagnetic 
dipole radiation. In a more general spin-down model, the relation between Cl and D. is deter- 
mined by a power law as 

Cl = -ka", (2.5) 

where A; is a constant and n is the breaking index. The integration of spin-down law of Eqn. 2.5 
yields 

where t is the time that passes until the angular velocity has changed from the initial value D-q to 
the current value Cl. This equation can be used to estimate a pulsar's age. With the reasonable 
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assumption 3> ^, the last term can be neglected and one obtains the characteristic age (t) 
as 

1 a _ 1 / _ 1 P 

where / and P are the frequency and period of rotation. 

Since for many pulsar n ^ 3 as expected for magnetic dipole breaking, n = 3 is assumed as 
a common definition for the characteristic ages, i.e. 

f P 

T=A = ^. (2.8) 
If IP 

Nevertheless, the breaking index can be correctly determined if Q. can be measured. It can be 
obtained by the differentiation of the spin-down law of Eqn. 2.5 and replacing k by the same 
equation as 

Similarly, one can determine the second deceleration parameter m through a second differenti- 
ation of Eqn. 2.5 as 



... n{ln-\)p m/3 ff 
/ = -ji = ^"^ = «(2«-l) = (2.10) 



according to which m is given through the third derivative of the pulsar frequency. Both deceler- 
ation parameters are therefore of high theoretical interest, since they reflect a pulsar's spin-down 
mechanism and deviations could indicate variations of the magnetic moment. 



2.2.4 Magnetosphere and Emitting Regions 

The magnetosphere is determined by the strength and the geometry of the magnetic field. A 
general estimate for the magnetic field {Bq) is given by 

5o = 3.3x 10i'^(P/')2G. (2.11) 

The magneto sphere's geometry is determined by the light cylinder — sometimes also called 
velocity-of-light cylinder. The light cylinder is a cylinder oriented parallel to the axis of rotation 
of the pulsar and its radius is given by the distance — c/Q. dX which a co-rotating particle 
would exceed the velocity of light (c). The magnetic field lines, which extend beyond this 
cylinder, are called open field lines, and the others closed. Inside the light cylinder, the pulsar 
contains a high energy plasma, as well as a relativistic stream of electrons and positrons — 
most of which co-rotate with the pulsar as they are confined to the field lines. Outside this 
cylinder the particle density is much lower, and the particles cannot maintain co-rotation due 
to the limit imposed by the velocity of light. Therefore, co-rotation is limited to particles on 
closed field lines inside the light cylinder, while particles and field lines at greater distances 
are wrapped in spirals around the pulsar. This geometry and the strong magnetic fields lead 
to two distinct regions in the magnetosphere where the radiation of the pulsar is generated: 
the polar cap and the outer gap. Both regions produce very high energy particles by different 
acceleration mechanisms and thus exhibit different radiation characteristics. Evidence for the 
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Figure 2.12: Magnetosphere of a pulsar showing the regions of the polar cap and outer gap which 
are presumably the regions of high particle acceleration and emission of radiation. (Figure taken 
from GiuUani [2006].) 



polar cap model is seen in beamwidth and polarization, while evidence for the outer gap model 
is seen in the high energy radiation from young pulsars. 

The polar cap region is indicated for one pole in Fig. 2.12. It is a small vacuum region near 
the magnetic poles in which charged particles can be accelerated. This region is believed to be 
the origin of highly polarized narrow radio beams. The radiation is coherent, and also circular 
polarization is often observed as the dominant component of polarization. 

The outer gap is located far out in the magnetosphere, close to the velocity of light cylinder 
as shown in Fig. 2.12. It is assumed to be responsible for a very similar beam pattern which 
can extend over many decades of the electromagnetic spectrum as e.g. observed for the Crab 
Pulsar. The beam direction and the timing of the pulses is determined by the geometry of 
the magnetosphere near the light cylinder. Electromagnetic radiation is generated by charged 
particles which have been accelerated by the high electric fields in the charge-depleted region 
of the outer gap. Moreover, this radiation can also be amplified by cascades of pair creation 
(Sturrock[1971]). 
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2.3 Model of the Pulsar Wind Nebula 

If one assumes that MSH 15—52 and PSR B 1509— 58 are associated, this system has to be 
understood in the framework of PWNs. Therefore it is important to understand the concept 
of PWNs for a correct interpretation of its observations. This section gives an introduction to 
PWNs. Further details about recent experimental results can be found in the review of Gaensler 
and Slane [2006]. Further theoretical discussions are presented by van der Swaluw [2001]. 

The term pulsar wind denotes the flow of particles, which is generated by the pulsar. It 
mainly consists of relativistic and ultrarelativistic electrons which are ejected from the pulsar 
surface, or which are produced in cascades of pair creation in the emitting regions of the pul- 
sar's magnetosphere (Sec. 2.2.4). Therefore, the term electrons equally refers to electrons and 
positrons, i.e. to the leptonic component, in the context of pulsar winds. Also, the importance 
of the contribution from VHE hadronic particles to pulsar winds is debated. The hadrons are 
mainly nuclei which are emitted from the pulsar's surface or provided by the SNR. 

2.3.1 Evolution 

Since a PWN is typically embedded in an SNR, its evolution is determined by the evolution of 
the SNR. At early times <~1 kyr) the SNR expands freely at velocities > 5 — 10 x 10^ km/s. 
The expansion of the relativistic pulsar wind occurs rapidly and symmetrically within the un- 
shocked ejecta. This early stage is described by the magneto-hydrodynamic (MHD) model of 
Kennel and Coroniti [1984], which was initially developed for the Crab Nebula. This model 
discussed in the next section is likely to apply to MSH 15—52 as well, which is a similar young 
PWN. The ages of the Crab Nebula and MSH 15—52 are 1.0 kyr and 1 .7 kyr respectively. 

The evolution of a PWN takes a turn when the SNR evolves into the Sedov-Taylor phase 
(f >~1 kyr). At this stage, the SNR also develops a reverse shock in addition to its forward 
shock. The reverse shock first moves outward behind the forward shock and eventually moves 
inward, leading to a compression of the PWN typically at time spans of several kyr. It is 
assumed that the Vela PWN is in this stage (Fig. 2.13). The compression of PWNs by reverse 
shock has been modelled by Blondin et al. [2001]. Fig. 2.14 shows that these simulations can 
reproduce this evolution. However, MSH 15—52 does not yet appear affected by a reverse 
shock. For comparison, the Vela pulsar has an estimated age of ~ 10 kyr, which is about an 
order of magnitude above the age of the Crab Nebula and MSH 15—52. 

-40" ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^H 




Figure 2.13: Parkes 2.4 GHz map of the Vela SNR 
(Duncan et al. [1996]) showing the embedded Vela 
PWN. The cross indicates the location of the associated 
pulsar B0833— 45, while the white arrow indicates its di- 
rection of motion (Dodson et al. [2003]). The fact that 
the pulsar is neither located at nor moving away from 
the PWN's center suggests that a reverse shock interac- 
tion has taken place. (Figure taken from Gaensler and 
Slane [2006].) 
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Figure 2.14: Simulated evolution of a PWN in the reverse shock of an SNR. The reverse shock 
compresses the PWN to a fraction of its earlier size. Since the expansion of the SNR is simulated 
in a nonuniform medium with a density gradient (in the vertical direction), this simulation is able to 
reproduce the displacement of the PWN from the center as observed for the Vela SNR in Fig. 2.13. 
(Figure taken from Blondin et al. [2001].) 

2.3.2 The Model by Kennel and Coroniti 

The evolution of a PWN at early times (? <~1 kyr), when the SNR is expanding freely at 
velocities grater than 5 — 10 x 10^ km/s into the ambient medium, is described in the model by 
Kennel and Coroniti [1984]. It is a self-consistent, spherically symmetric IVIHD model which 
was initially developed for the Crab Nebula. It describes the flow of relativistic plasma and 
the magnetic field from the pulsar to the boundary of the nebula within an SNR. Hadronic 
components are neglected. The model distinguishes six different concentric regions of different 
astrophysical properties. They are illustrated in Fig. 2.15 and characterized in the order of 
increasing radius (r) as follows: 

• Region I is contained in a small region within the light cylinder (r < r^), where the pulsar 
wind is created. It is indicated by the small spot at the center of Fig. 2.15 and 2.16. 

• Region II extends from the light cylinder to the standing shock front (r^ < r < r^) which 
forms where the highly relativistic pulsar wind expands supersonically (v > c/^/3) into 
the SNR ejecta. This region is under-luminous. The corresponding region in Fig. 2.16 
extends from the light cylinder to the bright ring with a radius of ~ 10". According 
to Weisskopf et al. [2000] and Hester et al. [2002], the ring, which is located between 
the pulsar and the torus, corresponds to the shock front where the cold relativistic wind 
converts into a more slowly moving synchrotron-emitting plasma. 

• Region III contains the synchrotron-emitting plasma, which extends from the standing 
shock front to the boundary of the PWN (r^ < r < r;v> ^^v ~ 2pc for the Crab Nebula.) It 
represents the main portion of the PWN and emits most of the radio, optical and X-ray 
synchrotron radiation. In Fig. 2.16 this region corresponds to the torus of the PWN be- 
yond the inner ring. 
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Figure 2.15: Sketch of the theoretical MHD Figure 2.16: Chandra X-ray data of the Crab 

model by Kennel and Coroniti [1984] for the Nebula showing features predicted by the model 

Crab Nebula. The model predicts different re- of Kennel and Coroniti [1984]. For example the 

gions of the nebula. (Figure taken from Kennel inner ring corresponds to the shock front. (Fig- 

and Coroniti [1984].) ure taken from Weisskopf et al. [2000].) 



• Region IV lies outside the PWN but still inside the SNR. For the Crab Nebular this re- 
gion is at a distance of ~ 5 pc from the pulsar. This region is outside the visible area of 
Fig. 2.16. 

• Region V is the region of material swept up by the blast wave of the supernova. 

• Region VI lies outside the SNR (r > 5 pc for the Crab Nebula). It only contains the 
interstellar medium. 

The experimental evidence from recent Chandra X-ray observations for the existence of 
these theoretically predicted regions is discussed in greater detail by Weisskopf et al. [2000] 
and Hester et al. [2002]. The agreement between Fig. 2.15 and Fig. 2.16 from Weisskopf et al. 
[2000] is remarkable. In particular the bright ring between pulsar and its torus is conspicuous. 

2.3.3 7-Ray Production in PWNs 

The model by Kennel and Coroniti [1984] can explain the observed radiation from PWNs of 
leptonic origin. As the regions differ in their astrophysical properties, they also differ in the 
production mechanism of electromagnetic radiation. IVlost important is the contribution of the 
inner three regions (I- III), which are illustrated in Fig. 2.17. Although the radiation of PWNs 
extends from radio to TeV y-ray energies, this section will mainly discuss the production of 
VHE 7 radiation. 

• In region I, i.e. inside or close to the light cylinder, the / emission is dominated by pulsed 
curvature, synchrotron and IC radiation. According to the polar cap model, electrons 
with an energy of ~ 10 TeV produce y-rays of ~ 10 GeV At 10 GeV a cutoff in the y-ray 
spectrum is expected, since the optical depth increases drastically with the 7 energy, and 
the radiation from the inner magnetosphere is heavily absorbed due to the strong magnetic 
fields. 
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2.3 Model of the Pulsar Wind Nebula 



• Region II is the regime of an ultrarelativistic, cold and under-luminous pulsar wind. The 
reduced luminosity is due to a reduction of curvature and synchrotron radiation. While 
it is obvious that curvature radiation decreases with decreasing curvature of the magnetic 
field lines, the decrease of synchrotron radiation can be explained by the fact that the 
wind and the magnetic field move together, since the field is frozen into the wind. On the 
other hand, an increase of IC radiation is unavoidable due to scattering of photons e.g. 
from the cosmic microwave background radiation (CMB), the interstellar radiation field, 
the emission from the magnetosphere and the thermal emission from the pulsar surface 
(~ 10^ K) (Bogovalov and Aharonian [2000], Aharonian and Bogovalov [2003]). The 
IC spectrum is primarily determined by the wind's Lorentz factor. Typical values of the 
latter range from 1 0^ to 1 0^ and result in IC y radiation from 1 GeV to 1 TeV. 

• In region III the shocked pulsar wind produces synchrotron and IC radiation. Since this 
region is much larger than the previous regions, it contributes the largest fraction of the 
observable emission from PWNs. Also, if a VHE hadronic component contributes to the 
/production through ;r^-decay, its main contribution would be expected from this region. 

Although the confirmation of 7 radiation of hadronic origin from decay appears more 
difficult, a contribution from hadrons is also discussed. Such considerations are interesting 
since they, in comparison to leptonic models, explain a larger fraction of the unexplained energy 
loss rate of pulsars which is derived from their spin-down luminosity (Horns et al. [2006]). 
While the leptonic wind observable by its IC radiation only transports a smaller fraction of 
the pulsar's spin-down luminosity, the acceleration of a nucleonic wind could absorb a larger 
fraction. The nucleonic interaction regions are inside the nebula and at the boundary to the 
surrounding interstellar medium. 

Radiation from a Pulsar-wind-nebula complex 




Figure 2.17: 7-ray production by electrons 
and positrons (e) from a pulsar wind. The 
schematic representation shows the different 
production mechanisms in the inner three re- 
gions of a PWN as described by the model 
by Kennel and Coroniti [1984]. R, O, X 
and 7 stand for radio, optical, X-ray and 
7-ray emission, respectively. CR, IC and 
Sy stand for curvature, inverse Compton and 
synchrotron radiation, respectively. The ori- 
entation of the magnetic field lines (B) is also 
indicated. (Figure taken from Aharonian and 
Bogovalov [2003].) 



IC 
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Interstellar medium 
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2.4 Radiation Mechanisms 

Various very high energy radiation mechanisms can explain the energy spectrum from PWNs 
observed. The important leptonic mechanisms are synchrotron, curvature and inverse Compton 
radiation. Leptons in this case again refer to electrons and positrons. The hadronic radiation 
mechanism is mainly the decay of K^, which are produced in interactions of protons and nuclei. 
Here the individual mechanisms are introduced and compared. 



2.4.1 Synchrotron Radiation 

Synchrotron radiation does not only dominate in many PWNs but also in many high energy 
astrophysical processes. Although in principle any charged high energy particle can emit syn- 
chrotron radiation in a magnetic field, the main contribution is observed from electrons and 
positrons. Synchrotron radiation is emitted when a charged particle with relativistic energy is 
forced to a spiral path by a magnetic field. A detailed discussion can be found in Ginzburg and 
Syrovatskii [1965], Blumenthal and Gould [1970] and [Longair, 1994, Chp. 18]. Here only a 
few important aspects are presented. 

The direction of the emission depends on the pitch angle 6, i.e. the angle between the 
velocity of the particle with respect to the field lines. The mean energy loss of an electron (or 
positron) is 

=2GTcUmag{iyfsm^e, (2.12) 

where c and v are the velocity of light and the electron respectively and 7 is the Lorentz factor 

8;r 

T 

is the Thomson cross section, is the classical electron radius, 

Umag = :r- (2.14) 
2jU0 

is the energy density of the magnetic field B and jUq is the permeability of free space. For an 
isotropic distribution of the pitch angles, the mean radiation loss of an electron is 

"(^) =l<'TcUmagQ'f = e.6xl0yB^cys-\ (2.15) 

The intensity spectrum emitted by a single electron can be written as 

= BsinGF — , (2.16) 



(7r = —r^ = 0.665 X lO-^V^ (2.13) 



dEg Sn'^eocnie V 



where e and trie are the charge and mass of an electron, 80 is the permitivity of free space, is 
the critical frequency for synchrotron radiation and F(x) is given by 



F{x)=x i^5/3 {z)dz ^ I (2. 17) 

(f)^x2exp(-x), ifx>l, 



18 



CHAPTER 2 MSH 15-52 and PSR B1509-58 



2.4 Radiation Mechanisms 



where ^5 p and T denote the modified Bessel and the Gamma function respectively. The critical 
frequency is given by 

Vc = ^{^)fvgsme, (2.18) 
where Vg is the gyrofrequency. Therefore, 

'dEe 



, , 0.29/1 Vc, (2.19) 

I sy 

where h is Planck's constant. The corresponding intensity distribution is shown in Fig. 2.18. 

The energy spectrum of synchrotron radiation for a sample of electrons, which has an energy 
distribution according to a power law, i.e. 

^=,c£-«, (2.20) 
dEe 

where E^ is the electron energy, jc is a normalization constant and a is the constant of the 
spectral index, results in a y-ray spectrum (<I>y) of the form 

d^y («+l) -(«-') 

— ^oc^5^£y^, (2.21) 

dEy ' 

where Eg is the energy of the photon and B is the magnetic field. So in this common case the 
relation between the photon index (F) of the 7-ray spectrum and the spectral index of parent 
electron distribution (a) is 

F=-^. (2.22) 
The cooling time, i.e. the time for a particle to lose all its energy, is calculated as 

Ee 



dt 



(2.23) 



The formalism of synchrotron radiation from protons with energy Ep is identical and a com- 
parison can be found in [Aharonian, 2004, Sec. 3.3.2]. Due to the larger proton mass (m^ = 
1836me), the radiation loss and therefore the cooling times from electrons to protons differs by 
orders of magnitudes. The ratio is given by 

/dEe\ ^ ^ 5 

^ = ^:^=f^V = 1.5xlO«. (2.24) 

dEp\ \me I 



dt , 

P 



Therefore synchrotron radiation from electrons dominates in most high energy astrophysical 
processes. 



2.4.2 Curvature Radiation 

Curvature radiation is similar to synchrotron radiation. It is also caused by acceleration of 
charged particles that pass through a magnetic field. In contrast to synchrotron radiation how- 
ever, where the acceleration due to gyration is perpendicular to the field lines, curvature radia- 
tion is associated with the acceleration parallel to, i.e. along, the field lines. Although curvature 
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radiation is usually exceeded by synchrotron radiation by orders of magnitude, curvature radi- 
ation becomes relevant for extremely strong and magnetic fields, such as those in the vicinity 
of pulsars. Since the radiation process is in many aspects similar to that of synchrotron emis- 
sion, the same equations apply, if the cyclotron radius r^y = jmv/eB is replaced by the radius of 
curvature p of the magnetic field lines. Then, similar to Eqn.2.12, one finds the radiation loss 
given by 

'iEi) =pLtr:)\ (2.25) 

dt eneop^ \cJ 

The intensity spectrum follows the distribution of synchrotron radiation given by Eqn. 2.16 if 
the critical frequency for curvature radiation (Vc) is used, i.e. 

^- F -] , (2.26) 



dEj STteop 

where 

v. = |^-lr'. (2.27) 

Therefore, Fig. 2.18 also represents the intensity spectrum of curvature radiation. 

A sample of electrons with a power law energy distribution given by Eqn. 2.20 results in a 
7-ray spectrum with a photon index of 

r=-^ (2.28) 
(cf. [Lyne and Graham-Smith, 1998, pg. 180]). 



2.4.3 Inverse Compton Radiation 



Inverse Compton (IC) radiation is produced in the scattering process of high energy particles 
with photons. The energy (Ey) of a photon after scattering with an electron at rest is given by 



Ey{Epii, 6) — EphP{Eph, 6), 



(2.29) 



where is the scattering angle, Eph is the initial energy of the photon and P{Eph, Q) is the ratio 
of the photon energy after and before the collision given by 



P{Eph,d) = 



1 



1 + -H(1-COS0) 



(2.30) 



where is the mass of an electron ([Longair, 1992, pg. 99]). The cross section (okn) is 
according to the Klein-Nishina formula 



^ = 1^ [piEph,e)-PiEph,e)hm^ie)+PiEph,e)'] , 



(2.31) 



where rg is the classical electron radius, c is the velocity of light. 

From these principles Blumenthal and Gould [1970] calculated the IC energy spectrum 
which is produced in the interactions of accelerated electrons with photons as 



d^y _ Inr^meC^ n{Eph)dEph 



dEy 



X 



(2.32) 
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where the constants p and q are given as 

mgc^ p(l-£e) 

/is the Lorentz factor and Eg is the initial electron energy. Fig. 2.19 shows the corresponding 
intensity distribution. 

Thomson scattering is obtained if Ey{Epii, 6) ft; Epfj, i.e. 

P{Eph,e)^l^^^{l-cose)<l. (2.34) 

In IC scattering the electron is not at rest but has the Lorentz factor of 7. Then the energy of the 
photon is ^ 2YEpfi (for v f=i c) in the rest frame of the electron. Therefore one distinguishes 

r ^yEph ^i^y^!^ jhe Thomson limit, 
if < ^ fn (2.35) 

4^ >l^Y>Tf- the Klein-Nishina limit. 

In the Thomson limit the maximum (Ej) and mean energy (Ej) of the scattered photons are 

Ey^4fEph (2.36) 

and 

1 T _ / Ep 



Ey fti -Ey Ri fEph = j Eph, (2.37) 

where Eg is the energy of the electrons. The Thomson limit is valid for many astrophysical 
processes. For example, for the CMB {Eph = 2.35 x 10""^ eV) the Thomson limit is fulfilled for 

^<i^-4^^i^ = «-5x>»' Ee = rn>.c^<250TeV. (2.38) 

In the Thomson limit the radiation loss is 

where Uph is the energy density of the photon field. The ratio of IC to synchrotron radiation is 
immediately given by the ratio of the photon to the magnetic energy density as 

* //c _ Uph ^2.40) 



dt 

I sy 

According to Ginzburg and Syrovatskii [1965], IC radiation in the Thomson limit of parent 
particles with an energy distribution according to a power law (Eqn. 2.20) passing through a 
monochromatic photon field has a photon index of 

r=^. (2.41) 
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Figure 2.18: Intensity spectrum of synchrotron 
and curvature radiation emitted by a charged 
particle of a fixed energy in a magnetic field of 
constant strength and curvatures. The spectrum 
is shown as a function of x = v/Vc, where Vc is 
the critical frequency. (Figure taken from [Lon- 
gair, 1994, pg. 247].) 




Frequency in units of t'/^'i^ 

Figure 2.19: Intensity spectrum of inverse 
Compton radiation emitted by a charged parti- 
cle. Vo is the frequency of the unscattered pho- 
ton. (Figure taken from [Longair, 1992, pg. 
104].) 



The radiation loss in the Klein-Nishina regime is 

According to Blumenthal and Gould [1970], IC radiation in the Klein-Nishina limit of parent 
particles with an energy distribution according to a power law (Eqn. 2.20) passing through a 
monochromatic photon field has a photon index of 

r=a+l. (2.43) 

2 4 

Calculations for IC radiation in the intermediate regime at relativistic energies (7~ can 
be found in Aharonian and Atoyan [1981]. 

Although high energy protons can also produce IC radiation, IC scattering of protons with 
the same energy as electrons is suppressed by many orders of magnitude as 

( M^ = 9x 10-14 (2.44) 
\mpj 

(cf. [Aharonian, 2004, Chp. 3.2]). Therefore, IC radiation from electrons dominates in high 
energy astrophysics. 



2.4.4 Hadronic 7 Radiation 

/radiation from hadronic interactions is mainly produced through the decay of secondary parti- 
cles from inelastic nucleon collisions. While the secondary particles are mainly pions, i.e. K^, 
K~ and with equal probability, only the K*^ mesons decay into two photons and contribute 
to the 7-ray spectrum. The mesons decay into muons, electrons and neutrinos. The ma- 
jority of nucleonic interactions are produced by highly accelerated protons which collide with 
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Figure 2.20: Total cross section for inelastic 
proton-proton collisions. The threshold energy 
^th is nipC^ + 2mjtC^{\ +m„/Amp), where the 
kinetic energy of the proton exceeds ~280 MeV, 
i.e. twice the mass of the tt". (Figure taken from 
Kelneretal. [2006].) 




Figure 2.21: Energy spectra of y-rays and other 
secondary products in proton-proton collisions. 
The proton energy spectrum obeys the power 
law of Eqn. 2.48. The electron and neutrino 
spectrum coincide. (Figure taken from Kelner 
et al. [2006].) 



ambient hydrogen, i.e. by proton-proton collision. Therefore, the inelastic part of the total 
proton-proton cross section (CineK^^p)) determines the hadronic y-ray spectrum. According to 
Kelner et al. [2006] C7inei(£^p) is approximated as 



Oinei(^/>) = (34.3 + 1.88L + 0.25L2) x 



mb, 



(2.45) 



where Ep is the energy of the proton, L = ln(£'p/l TeV) and is the threshold energy of 
the proton for the production of mesons. Since the kinetic energy of the proton has to 
exceed twice the rest mass of the pion (mj^ = 135 IVIeV) Eih = nipC^ + 2mjiC^{l +mjt/4mp) = 
1.22 X 10^3 TeV. Oinei{Ep) is shown in Fig. 2.20. 
The y-ray energy spectrum is then given as 



dEy 



qy{Ey 



dEj^^ 



(2.46) 



where Emin 
pions with 



Ey + m^c'^ /AEy, E-ii and qn{ETi) are the energy and the emissivity of secondary 



CUh 



-Cinel 



nipC^ + — ]Jn 



nipC 



(2.47) 



Here hh is the density of the ambient hydrogen, c is the speed of light, K is the mean fraction 
of kinetic energy of the proton transferred to 7 photons or the mesons per collision and Jp is 
the energy spectrum of the protons (Aharonian [2004]). 

A distinct feature of the hadronic 7-ray spectrum is a pronounced peak in the energy spec- 
trum at £■ = m;rC^/2 ~ 68 IMeV, which is independent of the energy distribution of the mesons 
and therefore of the proton. Fig. 2.21 shows the energy spectrum of the y-ray photons and the 
other secondary particles for protons with a spectral distribution 



pa 



exp 



Ep\ 



(2.48) 
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E, TeV E, TeV 

Figure 2.22: Spectral energy distribution corresponding to Fig. 2.21 (left) and for a harder proton 
spectrum with a spectral index a = 1 .5 (right). The latter spectrum is closer to the nucleonic emission 
model for MSH 15—52 by Bednarek and Bartosik [2003]. (Figure taken from Kelner et al. [2006].) 



where Ep is the energy of the proton, the spectral index a = 2, the cut-off energy Eq = 1000 TeV 
and /3 = 1. The peak in the y-ray spectrum at ~ 68]VIeV is apparent. Fig. 2.22 shows the 
corresponding SED for the same spectral index and for a different index of a = 1.5. Proton 
spectra with a harder index (a = 1.5) become important in the models of hadronic particle 
acceleration in PWNs by Bednarek and Bartosik [2003]. 

The average energy loss of a proton of about 50% at each collision is described by the 
coefficient of inelasticity (/ = 0.5). Taking this energy loss into account and assuming an ap- 
proximately constant ainei at high energies as Aharonian [2004], the y-ray spectrum reproduces 
the proton spectrum, and the photon index reads 

r ^ a. (2.49) 

Moreover, for a hydrogen density of no = n/ cm^ the cooling time is 

Tpp = ^ ^5.3 X 10^ yr (2.50) 

noOppfc n 



and 



dE„ \ E„ 

-jr) =-^ = mOppfcE. (2.51) 

/ pp ^pp 



2.4.5 Energy Spectra from PWNs 

The radiation mechanisms discussed above can provided important information for the under- 
standing of PWNs. Since they lead to different energy spectra from a common primary electron 
distribution, the radiation at different wavelength has to result in a consistent picture of the as- 
trophysical conditions at the source region. The different photon indices F for the same power 
law electron distribution with index a are summarized in Tbl. 2.3. 

Moreover, for a typical magnetic field strength of 10^^ G in a PWN the electron energy Eg 
for producing synchrotron radiation is approximately given by 

Ee = {10TcY)Bi^El^, (2.52) 
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were 5_5 = B/{10^^G) is a number for the transverse magnetic field strength (B) and Ey^^y is 
the number for the mean energy of the synchrotron radiation in units of keV (de Jager [2006a]). 
The corresponding value for IC radiation produced by scattering of CMB photons is approxi- 
mately given by 

£e=(18TeV)4eV' (2-53) 

where ExeV is the number for the mean energy of the IC radiation in units of TeV (de Jager 
[2006a]). From these two equations one immediately obtains the relation between the magnetic 
field strength, the synchrotron and the IC radiation as 

Ekev = 0.065:2£TeV- (2.54) 

It allows to infer the magnetic field strength if the synchrotron and IC component are both 
known. 

Also, one can explain spectral steepening with increasing distance form the center of ex- 
tended PWNs if lifetimes are considered. For example (cf. Aharonian et al. [2005a], de Jager 
[2006a]), the lifetime T(£'y) of VHE y-ray emitting electrons in a magnetic field is given by 

T{Ey) = (4.8kyr)5:2£Tei- (2-55) 
The corresponding lifetime t(£'x) for keV emitting electrons is 

t(£x) = (1.2kyr)5:|£;i. (2.56) 

In both cases the lifetime of electrons with higher energy is shorter. So with increasing dis- 
tance less high energy electrons and therefore a steeper photon index is expected as observed in 
Aharonian et al. [2006b]. 

Nevertheless, the modelling of astrophysical processes in a PWN often remains a difficult 
task, since many parameters are often not well constrained allowing for different explanations. 
However, simulations can often provide plausible solutions. 

Fig. 2.23 shows the energy spectrum of the Crab Nebula. It is an example of a PWN with 
a well determined energy spectrum over more than 20 decades. The synchrotron peak at keV 
energies and the IC peak at TeV energies are visible. The corresponding energies of the parent 
electrons producing the radiation are indicated by the labeled arrows. 

The energy spectrum of MSH 15—52 is less well constrained by measurements. However, 
different spectra have been predicted. A few are shown in Fig. 2.24 as calculated by Bednarek 
and Bartosik [2003]. These calculations show the contributions of the leptonic and also of 
nucleonic components for different densities of the ambient medium. The nucleonic energy 
spectra are similar to those used in the calculations shown in Fig. 2.25. They represent the equi- 
librium spectra of different nuclei after 1 kyr. The corresponding y-ray spectra are represented 
in Fig. 2.24 by the thin and thick dot-dashed curve. 
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Table 2.3: Relation between the photon index (F) and spectral index of a power law parent particle 
distribution (a) for different radiation mechanisms. 



synchrotron 


curvature 


IC radiation 


IC radiation 


7t^ decay 


radiation 


radiation 


Thomson limit 


Klein-Nishina limit 




r (a-l)/2 


(a-l)/3 


(a+l)/2 


a+l 


a 




Lo«(v/Hz) 



Figure 2.23: Non-thermal radiation ob- 
served from the Crab nebula. The solid and 
dashed curves represent the synchrotron 
and IC component radiation calculated for 
a spherical symmetric MHD model. The 
vertical arrows indicate the contribution to 
synchrotron radiation which is produced 
by parent electrons with energies of the 
corresponding label. (Figure taken from 
[Aharonian, 2004, pg. 47].) 
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Figure 2.24: Non-thermal radiation from 
MSH 15-52. The dashed and dot-dot- 
dot-dashed curves show the calculated syn- 
chrotron and IC spectra, respectively. The 
dotted and full curves show the same cal- 
culations with additional infrared photons. 
The nucleonic radiation from 71° decay for 
the density of the interstellar medium of 
0.3 cm~^ (thin dot-dashed curve) and for a 
high density of 300 cm^^ (thick dot-dashed 
curve) is also shown. (Figure taken from 
Bednarek and Bartosik [2003].) 
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Figure 2.25: Parent energy spectra of nu- 
clei of different mass numbers (A) in a 
PWN after Ikyr. Similar spectra have 
been used for the calculations in Fig. 2.24. 
The curves represent: A = 1 (dot-dot-dot- 
dashed curve), 2-10 (dotted), 11-20 (dot- 
dashed), 21-40 (dashed) and 41-55 (full). 
The corresponding 7-ray spectra are shown 
in Fig. 2.24 for different densities (thin 
and thick dot-dashed curve). (Figure taken 
from Bednarek and Bartosik [2003].) 
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7-Ray Astronomy with Imaging 
Atmospheric Cherenkov Telescopes 

Astroparticle physicists and astronomers have long been interested in observing cosmic VHE 
y radiation. After many years of research, imaging atmospheric Cherenkov telescopes (lACTs) 
were developed and established as useful instruments in VHE 7 astronomy. lACTs can detect 
cosmic / radiation through its interaction with the earth's atmosphere. Since the lACT technique 
allows for building systems with large effective areas, it is possible to detect a significant amount 
of cosmic / radiation in the energy range from about lOOGeV to lOOTeV. At this energy range 
even the strongest sources have a flux of less than 10^*^cm^s ^ LACTs are currently the most 
sensitive instruments for 7-ray astronomy at these energies. 



3.1 Air Showers 

To understand the lACT technique it is useful to know about the physics of air showers, which 
develop in the atmosphere, and about the Cherenkov light which is emitted. One can distinguish 
between two types of showers, electromagnetic and hadronic showers. 

3.1.1 Electromagnetic Showers 

Electromagnetic showers are electron-photon cascades which are initiated when a photon of 
high energy enters the atmosphere. The interaction with the molecules of the atmosphere leads 
to pair production. In turn, the produced electron-positron pairs emit photons via bremsstrah- 
lung. The sequence of pair production and bremsstrahlung results in a cascade with an ex- 
ponential increase of particles. The maximum of particles is reached when their energy has 
reduced to critical energy (Ec). At the particles' energy is not sufficient to sustain the pair 
production, and the remaining energy is finally dissipated by ionization. A full shower cascade 
develops within ^-^50 micro seconds. The frequency of the interactions is determined by the 
radiation length Xq. Xq is defined as the mean distance the particles travel when they lose all 
but exp( — 1) of their energy. Similarly, the interaction length (Aq) is defined as the distance a 
particle traverses until the probability is exp( — 1) that no interaction will occur. The interaction 
length for pair production is about | of Xq for bremsstrahlung. After crossing the distance X a 
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particle's energy is given by 

E{X) = Eooxp(-^], (3.1) 



^0. 

where Eq is the initial energy. Bremsstrahlung and pair production mainly happen in interactions 
with the nuclei of the atmosphere, since the probability is proportional to the square of the 
atomic number. 

In a simple model by Heitler [1954] (Fig. 3.1), the differences between the radiation and 
the interaction length are neglected and it is assumed that bremsstrahlung and pair production 
occur with the same frequency. Then, the distance R, after which on average half the particles 
interact, is defined through 



Jo J 2 

^R = Xoln2. (3.3) 

The total number of particles after n steps of interaction isN = 2". Assuming that the particles 
split their energy at each interaction, the critical energy is reached after the maximal number of 
{nmax) interactions. Therefore, 

En 

Ec = ^, (3.4) 
which yields the relation between Umax Ec as 

_ HE^/Ec) 

''■max — , • W'-'/ 

In 2 

Moreover, one obtains the relation for the number of particles at the shower maximum Nmax by 

En 

Nntax = 2"'^ = -^, (3.6) 

with a composition of ^N^ax photons and ^N^ax electrons and positrons. Also, the depth of the 
shower maximum ?njax the atmosphere is given by 

^max — ^max^Oi (3-7) 

where is the radiation length measured in matter per cm^. = pXo, where p is the density. 

With the values of Ec and for air, one can calculate these shower parameters for the 
different energies of the primary particle. A few values are given in Tbl. 3.1. In comparison with 
the more accurate Monte Carlo simulation of the longitudinal shower development in Fig. 3.2, 
these values already provide a good estimate. A more accurate description is given by the 
Nishimura-Kamata-Greisen (NKG) formula (Greisen [1965]). 

The height of the shower maximum depends on the atmospheric depth which is determined 
through the atmospheric density profile 

p(/i) = 1.3 X lO-Vm-^exp (^^^ . (3.8) 



i^For air = 87MeV (Wigmans [2000]). 
2)For air = 37 gcm-^ (Yao et al. [2006]). 
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3.1 Air Showers 




height [m] (for vertical showers) 




Figure 3.2: Average longitudinal 
profile of electromagnetic air show- 
ers with energies for the primary 
y-ray interaction between 0. 1 TeV 
and 1 PeV. N^. is the number of 
electron-positron pairs with energy 
above 5MeV as determined from 
Monte Carlo simulations. The rela- 
tion between the shower depth and 
shower height is given by the top and 
bottom scales. (Figure taken from 
Bemlohr [1996].) 



atmospheric depth [g/cm**2] 



Table 3.1: Parameters of electromagnetic air showers, rimax^ ^max and are the number of interac- 
tions, the number of particles and the depth at the shower maximum, respectively, tmax is obtained by 
a simple estimate using the Heitler model. The corresponding height {h) is determined from Fig. 3.2. 



Eq [TeV] 


f^max 


A^max[10=] 


tmax [Xq] 


h [km] 


0.1 


10 


0.01 


370 


10 


1 


13 


0.1 


500 


6 


10 


17 


1 


620 


4 



The relation between height and atmospheric depth is represented by the top and bottom scales 
of Fig. 3.2. 

The lateral distribution of a 7 shower is determined by Moliere scattering (Bethe [1953]), 
which describes the multiple Coulomb scattering of electrons and positrons in the atmosphere. 

The characteristic parameter is the Moliere radius. The spread caused by bremsstrahlung is 
negligible since bremsstrahlung is emitted in a cone in forward direction with an angle ~ 
= i. At sea level, 90% of the shower energy is deposited in a cylinder around the shower 
axis with a radius of 80 m. Fig. 3.5 and Fig. 3.6 show the lateral profile of a simulated shower 
of a photon of 300 GeV and the corresponding Cherenkov light at the ground, respectively. 
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3.1.2 Hadronic Showers 

Hadronic showers are initiated by nuclei of cosmic radiation (A^cosmic) which penetrate the at- 
mosphere. The cosmic radiation consists mainly of protons (87%), a particles (12%) and a 
small fraction of heavier atomic nuclei. Electrons, 7-rays and high energetic neutrinos only 
constitute a minor fraction. When entering the atmosphere the nuclei produce spallation frag- 
ments and new particles (X) in inelastic collisions with the nuclei of the atmosphere (A^atm.)- 
The new particles are mainly pions in the ratios ^K^, ^7t^, j^"- The have a lifetime of 
8.4 X 10^^^ s and decay into two photons which can initiate electromagnetic showers. The n"^ 
have a longer lifetime of 2.6 x 10^ s and can produce other particles, mainly pions, in inelastic 
collisions. The charged pions can also decay into muons and neutrinos. The muons decay into 
neutrinos and electrons, which can start electromagnetic showers. Therefore, a hadronic shower 
also has an electromagnetic shower component. Also the primary particles and the spallation 
fragments can form hadronic subshowers. So the main reactions are 



A'cosmic + A'atm. ^ hadrons+X (3.9) 

Jt^ r+r (3.10) 

7t- ^ fi + v^ (3.11) 

7t+ ^ ju+v^ (3.12) 

fi- e- + v, + Vju (3.13) 

H+ e+ + Ve + v^i. (3.14) 



The development of a hadronic shower is sketched in Fig. 3.3. Fig. 3.4 shows the distribution of 
high energy particles in the atmosphere. It also represents the composition of hadronic shower 
cascades which are the most frequent in the atmosphere. 



f 

I 




Figure 3.3: Model of a hadronic air shower. 
Characteristic is the hadronic component of pi- 
ons and muons. Electromagnetic subshowers 
are also part of a hadronic shower. (Figure taken 
from [Longair, 1992, pg. 149].) 




2000 4000 6000 9000 10000 

Depth through atmosphere/kg m-=^ 



Figure 3.4: The vertical fluxes of different com- 
ponents of particles in the atmosphere as pre- 
sented by A. M. Hillas (1972, Cosmic rays. Page 
50, Oxford: Pergamon Press.) (Figure taken 
from [Longair, 1992, pg. 150].) 
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3.2 Cherenkov Light 



The mean free path of the hadrons is given by 

A/ = — , (3.15) 
nOi 

where n and cr, are the particle density and the target particles' cross section in the medium. The 
atmospheric depth for hadrons corresponds to about 12 x A,. Since A, is about twice the radiation 
length of electromagnetic showers, hadronic showers penetrate deeper into the atmosphere. 

Due to the transversal momentum of the secondary particles, the hadronic showers have a 
larger lateral spread than electromagnetic showers, and are more irregular. Their lateral distri- 
bution can be used for the separation from electromagnetic showers. Fig. 3.5 and Fig. 3.6 show 
the longitudinal profile of a simulated proton shower of 1 TeV and the corresponding Cherenkov 
light at the ground, respectively. 



3.2 Cherenkov Light 

Air showers can be detected by their Cherenkov light. Cherenkov light is emitted from charged 
particles which travel faster through a medium than the speed of light in that medium. For a 
medium with an index of refraction n, the condition for the velocity v to produce Cherenkov 
light is 

v>-, ^ j8>-, (3.16) 
n n 

where c denotes the speed of light. This condition implies an energy threshold (Eth) which a 
particle has to exceed before it can emit Cherenkov light, namely 

Ethim) = Ythmc^ = ^ moc^ = ^ moc^. (3.17) 

Therefore, particles with a low mass, such as electrons, dominate Cherenkov emission. The 
threshold energies for electrons and muons are 21 MeV and 4.3 GeV respectively. 
Cherenkov radiation is emitted at an angle of 

1 

0c = acrcos— (3.18) 
np 

relative to the direction of the particle's velocity. For air with the index of refraction hq ~ 
1 + 3 X 10""^ and the condition of Eqn. 3.16, the maximal opening angle of 1.4° is obtained by 
j8 = 1, i.e. v — c. 

In air, Cherenkov radiation is emitted at wavelengths A between 400 nm and 700 nm. About 
30 photons per meter are produced by a single charged particle. Although it takes about 50 
microseconds for a shower cascade to develop, the front of Cherenkov light is only visible 
within ~10ns at the ground, since the cascade develops nearly along the light pass. At each 
point within the Cherenkov cone the light is only visible for 5 ns. Within 100 m from the shower 
axis, the light front reaches the ground with about 100 photons per m^. Therefore, lACTs 
require cameras with high sensitivity and short exposure times. 
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Figure 3.6: Monte Carlo simulations of Cherenkov light at the ground emitted from a / shower (left) 
and proton shower (right). The 7 shower shows the typical Cherenkov light pool with a radius of 
~120m. The proton shower is rather irregular. (Figure taken from Bernlohr [2006]). 
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Fig. 3.6 (left) shows the distribution of Cherenkov light at the ground, simulated for a 
7 shower of 300 GeV. The light is rather homogeneously concentrated within a radius of ~ 1 20 m 
from the shower axis and attenuates outside. This radius is similar for 7 showers of different 
zenith angles and energies. Fig. 3.6 (right) shows the same simulation for a proton shower of 
1 TeV. The light is inhomogeneously distributed and has several intensity maxima. 

3.3 Imaging Atmospheric Cherenkov Technique 

With lACTs it is possible to detect and identify the Cherenkov light from 7-ray showers in the 
atmosphere and thus to determine the energy and arrival direction of the 7-rays. Detection is 
possible with an optical system which consists of a mirror dish and a camera close to the focal 
plane. If a shower develops close to the telescope, the shower profile is reflected to the camera 
and recorded. "Close" means within the radius of the Cherenkov light pool of ~ 120 m at the 
ground. Fig. 3.7 illustrates the mapping of an air shower by the optical system of a telescope 
through geometric optics. 

Since the atmosphere is an essential part of the detection, lACTs are sensitive to the atmo- 
spheric conditions during observations. On the other hand, the collection areas do not have 
theoretical limits and they scale linearly with the mirror surface and the number of telescopes. 

lACTs can be combined into arrays for observation in a stereoscopic mode, meaning that 
a shower is recorded by at least two telescopes from different viewing angles. The advantages 
are e.g. a higher accuracy in the shower reconstruction and an improved rejection rate of back- 
ground events. 




Figure 3.7: Illustration of the imaging atmospheric Cherenkov technique. An air shower produces 
a light cone of Cherenkov light close to the center of a telescope array (left). The shower image is 
recorded by the telescope. The geometric optics of the mapping to the camera plane and the camera 
image are illustrated (right). (Illustration taken from [Schlenker, 2005, pg. 39].) 
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The separation of 7 showers from the majority of hadronic showers, the shower reconstruc- 
tion and data analysis are accomplished with Monte Carlo simulations of the air showers and 
the system's response. 

3.4 Overview of Current lACTs 

lACTs have been in use for about two decades now. The Whipple lACT is considered the 
first lACT which was widely recognized among y-ray astronomers. The Whipple collaboration 
established the lACT technique and detected TeV 7 radiation from the Crab Nebula in 1989 
(Weekes et al. [1989]). In subsequent years new lACTs were developed, among them Cat, 
HEGRA and Cangaroo. These are considered as lACTs of the first generation. They detected 
new TeV 7-ray sources and improved their technique. In 2002 H.E.S.S. came into operation 
which is considered as the first lACT of the second generation, for it is based on the same tech- 
nique but has a greatly increased sensitivity and precision. H.E.S.S. is currently one of the most 
sensitive lACTs and has more than doubled the number of known TeV 7-ray sources by 2006. 
Other similar sensitive lACTs of the second generation are Magic, Veritas and Cangaroo III. 

Geographic location is important for lACTs, since it determines the observable sky regions. 
Systems located in the southern hemisphere like H.E.S.S. have a direct view of the galactic 
center and the galactic plane where the majority of galactic TeV 7-ray sources are located. 
The map in Fig. 3.8 shows the second generation lACTs worldwide. The even distribution of 
telescopes between the northern and southern hemisphere is advantageous, since it provides full 
coverage of the TeV 7-ray sky. 




Figure 3.8: Imaging atmospheric Cherenkov telescopes worldwide. The even distribution of tele- 
scopes between the northern and southern hemisphere provides coverage of the full TeV 7-ray sky. 
Telescopes located in the southern hemisphere have a direct view of the galactic center and the 
galactic plane. (Map taken from Punch [2005].) 
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Figure 4.1: The H.E.S.S. array on the Khomas Highlands in Namibia. (Photo from Eifert [2005].) 

H.E.S.S. is the name of the new imaging atmospheric Cherenkov observatory located in 
Namibia (Fig. 4.1). It is an acronym for High Energy Stereoscopic System. H.E.S.S. was 
founded through the international collaboration of around 100 scientists from about 20 Euro- 
pean and African institutes under the leadership of the Max-Planck Institute for Nuclear Physics 
in Heidelberg, Germany. The name was also chosen in honor of the Austrian physicist Victor 
Franz Hess, who laid the foundations of modern astroparticle physics by his discovery of cos- 
mic rays in 1912 and who, as a result, was awarded a Nobel Prize in 1936. H.E.S.S. is an lACT 
of the second generation which succeeds the lACT experiment HEGRA. It has been developed 
as well as is operated partially by the same people. H.E.S.S. is very sensitive in an energy range 
from 0.2 to 50TeV. It is able to detect a y-ray point source that has a flux of 2.0 x lO^^cm^s \ 
corresponding to only 1% of the flux from the Crab Nebula with a significance of 5 C7 in about 
25 hours or a source of similar strength within 30 seconds (Aharonian et al. [2006a]). One of 
the basic concepts of H.E.S.S. is the technique of stereoscopy, as is reflected in its name. Stere- 
oscopy provides improved shower reconstruction and increased rejection rates for background 
events. In its first phase (H.E.S.S. I), the array consists of four identical Cherenkov telescopes 
(CT 1, CT2, CT 3, CT4). In the second phase (H.E.S.S. II), the array will be supplemented by 
an additional telescope located in the center of the array which will have a larger mirror surface 
and increased sensitivity. H.E.S.S. II is currently under development. H.E.S.S. I started obser- 
vation in June 2002 using its first telescope, and in the meantime the full telescope array was 
gradually completed. Since early 2004 observations have been made using the array of all four 

''The Crab Nebula is the standard candle in VHE 7-ray astronomy. 
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telescopes in stereoscopic mode. In 2006, H.E.S.S. had already confirmed most of the approxi- 
mately 10 TeV y-ray sources known before and had discovered about 20 new ones. Within the 
first two years of its operation H.E.S.S. exceeded most scientists' expectations (cf. Hofmann 
[2005] and Aharonian et al. [2005c]). 



The H.E.S.S. site is located in the Khomas Highland of Namibia (Fig. 4.2). The geographic 
location of the center of the telescope array is 16°30'00.8"E, 23°16'18.4"S at 1800m asl. 
There were several reasons why this location was selected. First, the dry climate of the Khomas 
Highlands allows for observations to be made throughout the year, with a total observation time 
of approximately 1600 hours per year at good and stable atmospheric conditions. The many 
km^ of sparsely populated area surrounding the H.E.S.S. site provide a minimum of night sky 
background also. Yet the capital of Namibia, Windhoek, lies at a distance of 100 km northwest 
from the H.E.S.S. site and thus provides the necessary infrastructure to maintain the observa- 
tory at reasonable costs. Finally, H.E.S.S.' location in the southern hemisphere permits the 
observation of the galactic plane and the galactic center at high zenith angles, which is when 
the telescopes are most sensitive. The galactic plane is particularly important for observation 
because it hosts most galactic y-ray sources and the super massive black hole Sgr A*. Besides 
the four Cherenkov telescopes, the site also contains the optical robotic telescope ROTSE, a 
control building with a workshop, a generator house and a residence building at a distance of 
1 km from the observatory buildings, separated from them by a hill. 



The four telescopes are placed in the corners of a square, which has a length of 120 m and its 
diagonals oriented in north-south and east- west directions. Each telescope is made with a rigid 
steel structure and has a total weight of 50 tons. Each consists of a reflector dish that is 13 m in 
diameter and a camera mounted in the focal plane of the telescope at a distance of 15 m from 
the dish. The dish is mounted at a height of 13 m above the ground to the support structure. The 



4.1 The Site 




Figure 4.2: Map of Namibia showing the 
location of the H.E.S.S. site. (Map taken 
from Schlenker [2005, pg. 40].) 



4.2 The Telescopes 



36 



CHAPTER 4 The H.E.S.S. Experiment 



4.2 The Telescopes 




Figure 4.3: H.E.S.S. telescope (CT 1) in the parking position. During the daytime the camera is 
protected in the shelter on the right. The scale is demonstrated by the person in the front. (Photo 
provided by Frank Breitling, 2004.) 



telescope can reach a maximum height of 28 m for observations at the zenith. Fig. 4.3 shows 
CT 1 in the parking position. During the daytime the telescopes are parked and the cameras are 
protected in their shelters against light, heat and dust. The scale of the telescope is demonstrated 
by the person standing in front on the structure. 



4.2.1 Davies-Cotton Design 

Each dish has a total reflector area of 107 m^. If shadowing by the camera support structure 
is taken into account the effective reflector area is reduced to about 95 m^. The reflector is 
composed of 380 circular mirror facets, each with a diameter of 60 cm and a reflectivity of 
80%. Each mirror is mounted on a support unit containing two actuators which allow individual 
alignment. The reflector follows a Davies-Cotton design (Davies and Cotton [1957]), which 
means that the dish and its mirror facets are spherical and have a focal length that is identical to 
the focal length of the dish. An advantage of the Davies-Cotton design is the cost efficiency of 
its manufacture. For in comparison to other designs, all of the mirrors in this one are identical. 
Although the Davies-Cotton design suffers from spherical aberration, this is not critical for 
H.E.S.S. The reason for this is that the residual point spread function of the reflector dish after 
the alignment of the individual mirrors is well contained in a camera pixel with a size of 0.16°. 
Also, the time dispersion due to aberration of approximately 5 ns is not critical for the readout 
of the camera image. Detailed information about the telescope mirror, its alignment and optical 
characteristics are given in Bemloehr et al. [2003] and Cornils et al. [2003]. 
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4.2.2 Pointing Accuracy 

Each telescope has an ah-az tracking with a slew speed of 100° min~^ The tracking position is 
measured by shaft encoders with a digital step size of 10" and is maintained with an accuracy 
of 30". Certain conditions can negatively affect the actual pointing position of the telescope's 
structure, mainly the camera support structure, which can bend due to its weight. But also 
wind pressure or dirt on the rail of the tracking system can eventually contribute to a loss of 
pointing accuracy of a few arc seconds. To monitor the deviations, each telescope is equipped 
with two CCD cameras: a Sky CCD and a Lid CCD. The Sky CCD is located in the right side 
of the dish and can monitor the telescopes field of view (FOV). The Lid CCD is located in the 
center of the dish and can monitor stars that are reflected by the dish onto the closed camera lid. 
From simultaneous observations by these two cameras, a pointing model has been developed 
(Gillessen [2004]) which describes the actual pointing position as a function of the tracking 
position. It is used to apply corrections off-line during the data analysis and is able to limit the 
systematic pointing error to 20". 

4.3 The Camera 

The H.E.S.S. cameras are described in great detail by Vincent et al. [2003]. A camera consists 
of 60 drawers which contain a total of 960 pixels. Each pixel consists of a photo-multiplier tube 
(PMT) with a quantum efficiency of 20-30% in the wavelength range of 300-700 nm. The front 
of the PMTs is equipped with a layer of hexagonal Winston cones in a honeycomb arrangement 
which reduces the light insensitive area between neighboring pixels to about 5%. The FOV 
of each pixel is 0.16° and contributes to the total FOV of 5° of the camera. Each PMT is 
calibrated to respond with an amplification of 2 x 10^ electrons for each collected photo-electron 
(p.e.). The PMT signal is fed into three different channels: the low gain, the high gain and the 
trigger channel. The low and high gain channels provide a linear response from 1 to 1600 p.e. 
Their signal is stored in an Analogue Ring Sampler developed by the ANTARES experiment. 
It samples the signal at 1 GHz over time windows of 16 ns. If an event has been triggered, 
the corresponding buffer is digitized and sent to the central data acquisition system. It takes 
about 610 jUs until the data is transferred and new data can be recorded. This is the dead-time 
of the camera. The resulting upper limit for the camera's acquisition rate is 1.6 kHz. The 
core camera's electronics are located in a crate behind the layer of PMTs. It contains, among 
others things, the sockets for the drawers, the readout and trigger cards, a central processor unit, 
the bus systems, a lOOMbits/s network interface, four power supplies, 16 temperature sensors 
and about 80 computer controlled fans. In addition, each camera is equipped with a global 
positioning system providing event times with an accuracy of /is. The electronics constitute 
the camera's data acquisition system. It is controlled by a Linux operating system written in 
programming language C. The camera's electronics and PMTs are housed in a container that is 
2mx2mxl.6m, with a total weight of about 900 kg and a lid in the front and in the back. Its 
total power consumption is about 5 kW. 

While the individual pixels constitute a first level trigger, the camera trigger system as a 
whole constitutes the second level trigger. It consists of 38 overlapping trigger sectors, with each 
containing 64 pixels. The typical trigger condition requires four neighboring pixels to exceed 
a threshold 5 p.e. within a time window of 2 ns. It takes about 70 ns to build a trigger signal 
which is fast enough to read out the data from the Analogue Ring Sampler. The camera trigger 
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Figure 4.4: Views of the second H.E.S.S. camera. The rear view (left) shows the crate with the four 
power supplies, the bus system and the network interface. The front view (middle) shows the 960 
Winston cones. When the camera has reached observation position, the lid opens for observations 
(right). (Figures taken from Vincent et al. [2003].) 



significantly reduces the number of background events. Depending on the trigger configuration 
and the zenith angle of the observations, the typical camera trigger rate is about (200 ± 50) Hz. 

4.4 The Central Trigger System 

In addition to the individual camera triggers, H.E.S.S. has a central trigger system (CTS) which 
constitutes a third level trigger. A detailed description of this system is given by Funk et al. 
[2004]. The CTS is designed to identify stereoscopic events by coincidence of individual tele- 
scope triggers. The standard central trigger condition requires a minimum of two telescope 
triggers. Higher trigger multiplicities provide events of higher reconstruction quality but at a 
cost of a reduced sensitivity for events of lower energy. The stereoscopic trigger condition ef- 
fectively reduces the background, e.g. muon events, which trigger the cameras. In addition, the 
CTS reduces the dead-time of the individual cameras, sending reset signals to triggered cam- 
eras to stop the readout process if no coincidence with other telescope triggers is observed. The 
CTS also measures the system's dead-time as well as assigns unique numbers to events which 
allows for individual camera images to be combined as a single event. Depending on the zenith 
angle of the observation, the typical camera trigger rate is about (300±50) Hz. All of the CTS's 
hardware is located in a crate (Fig. 4.7) in the control building. 

4.5 Atmospheric Monitoring 

H.E.S.S. has instruments for monitoring the atmosphere and for providing specific information 
about atmospheric conditions during its observations. A detailed description of this is given by 
Aye et al. [2003]. Atmospheric data is displayed in the control room and informs the observers 
about the applicable atmospheric conditions. It is also sent to the central data acquisition system 
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and recorded for off-line data analysis. The atmospheric data provides important information 
for data quality selection. The different monitoring devices involved are briefly described in the 
following subsections. 

4.5.1 Radiometers 

H.E.S.S. has five radiometers — one on every telescope and one scanning radiometer in the 
center of the array. The telescope radiometers measure infrared radiation from the sky in the 
FOV in a transmission window between 8 to 14/jm and calculate the temperature of the atmo- 
sphere through comparison with black body radiation. Since clouds reflect ambient light, their 
spectrum differs significantly from the spectrum of the clear sky and thus clouds can easily be 
detected. The scanning radiometer works the same way but scans the whole sky for the presence 
of clouds and approaching weather fronts. 

4.5.2 Ceilometer 

The ceilometer (Fig. 4.5) consists of a LIDAR (Light Detection and Ranging) system which 
emits short laser pulses (Brown et al. [2005]). It is located next to the scanning radiometer. By 
measuring the amount of backscattered light and time of flight, it can provide a detailed density 
profile of aerosol in the atmosphere and detect layers of clouds up to 7.5 km. The correlation 
between the amount of aerosol in the atmosphere and the trigger rate has been investigated by 
Le Gallou et al. [2005]. The aerosol absorption in the atmosphere affects all events, reducing 
the measured light intensity and thus shifting the energy scale. 

4.5.3 Weather Station 

Fig. 4.5 shows the weather station which is located in the center of the array. It consists of a ther- 
mometer, hygrometer, barometer, anemometer and pluviometer. Data from these instruments is 
continuously monitored, displayed in the control room and also recorded. 




Figure 4.5: Instruments at the H.E.S.S site for atmospheric monitoring: weather station (left) and 
ceilometer (right). (Photos taken from Schlenker [2005, pg. 46].) 
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4.6 The Central Data Acquisition System 



4.6 The Central Data Acquisition System 




Figure 4.6: Interior of the control room. The dis- 
plays provide the observers with the latest mon- 
itoring information. Terminals provide access to 
the central data acquisition system and control to 
the H.E.S.S. array. (Photo provided by Frank Bre- 
ithng, 2004.) 



The central data acquisition system (DAQ) is described in Borgmeier et al. [2001] and Borgmei- 
er et al. [2003]. It connects the individual H.E.S.S. components and manages the storage of data. 
It also provides the necessary interface for the control of the system. The DAQ is controlled 
through the DAQ front end in the control room (Fig. 4.6). Displays in the control room provide 
the observers with important monitoring information for observations such as trigger rates, cam- 
era images, PMT currents, significances of y-signals, the system load, weather and atmospheric 
information. Three computers provide the graphical user interface to schedule, configure and 
start observations as well as to stop them. The observations are typically scheduled in runs of 
28 minutes. 

The DAQ consists of object-oriented software in the computer language C-i-i-, which is based 
on the software packages ROOT (Brun et al. [2006]) and omniORB (Lo and Pope [1998]). 
During data taking, data is received from the individual components, then combined to events 
and stored in the ROOT file format. 

The software is run on a Linux computer farm consisting of about 20 Intel 386 compatible 
PCs, a gigabit network and two RAID arrays with a capacity of several terabyte each (Fig. 4.7). 
The hardware is located in the control building. New data is written to gigabyte tapes and 
shipped to European computing centers where it is calibrated and prepared for data analysis. 



Figure 4.7: The Linux PC 
farm in the control building 
is the core component of the 
central data acquisition sys- 
tem. The crate which houses 
the central trigger hardware 
is partially visible on the left. 
(Photo taken by Frank Bre- 
itling, 2004.) 
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After raw data has been recorded, it can be analyzed and scientific information can be extracted. 
This process of data analysis requires a good knowledge of the hardware, software and methods 
that can be applied in order to extract the desired information. The H.E.S.S. standard analysis 
has been developed to accomplish this task. It uses the well-established methods of imaging 
atmospheric Cherenkov astronomy which can reliably extract results from raw data in a com- 
putationally efficient way. The accuracy of the results has been analyzed in detailed Monte 
Carlo studies, through alternative analysis methods (Aharonian et al. [2006c]) and through di- 
rect comparison with results from other experiments (Aharonian et al. [2006a]). 

Several computational steps are necessary before physical quantities are obtained. Many 
of them require additional data from Monte Carlo simulations. The H.E.S.S. standard analysis 
(H.E.S.S. collaboration [2001]) performs these steps. It consists of a set of software packages 
based on the programming language C++ and the data analysis framework ROOT (Brun et al. 
[2006]). ROOT provides the statistical tools, which are used in the data analysis, e.g. for the 
fitting of functions. In the data analysis, some results rely on others so that all steps have to 
be performed in a certain order. This order is illustrated by the flowchart in App. A. Another 
discussion of the methods and the accuracy of the H.E.S.S. standard analysis can be found in 
Aharonian et al. [2006a]. 

The first part of the analysis consists of reconstructing y-ray showers. In a subsequent part 
statistical methods are applied, sky maps are generated and the energy spectrum is determined. 
In this chapter, the individual steps of data analysis are discussed and they are verified through 
application to H.E.S.S. data from the Crab Nebula. Sky coordinates are given in right ascension 
(RA) and declination (Dec) of the J2000 equatorial coordinate system unless stated otherwise. 



5.1 Monte Carlo Simulations 

Since many results of data analysis rely on data from Monte Carlo simulations, the later is 
essential for the analysis. For example, Monte Carlo simulations are required for calibration, 
background reduction, shower reconstruction and determinations of effective areas, i.e. deter- 
mination of the energy spectrum. Monte Carlo data for the standard analysis is simulated in two 
steps: first a shower simulation and then a detector simulation. 



43 



5. 1 Monte Carlo Simulations 



CHAPTER 5 The H.E.S.S. Standard Analysis 



5.1.1 Shower Simulation 

Shower simulations are produced with CORSIKA, a program for Cosmic Ray Simulations for 
Kascade (Heck et al. [1998]). CORSIKA can simulate atmospheric shower cascades and the 
corresponding Cherenkov light emission for various particles such as photons, protons and lep- 
tons. The simulations include many details, such as the effect of the local magnetic field of the 
earth on the charged particles in the shower cascade, in order to provide very realistic results. 
Also, the different atmospheric transmission profiles generated with MODTRAN (Anderson 
et al. [1996]) can be included in the simulation. For H.E.S.S simulations, two different atmo- 
spheric models are used which both describe the atmosphere at the H.E.S.S. site sufficiently 
well enough, as confirmed by trigger studies (Funk et al. [2004]). The desert model reproduces 
a clear atmosphere with little haze and a boundary layer starting at 1800 m above sea level. 
The maritime model reproduces a more humid atmosphere with a boundary layer starting at sea 
level. The desert model is generally preferred in most analyses. 

5.1.2 Detector Simulation 

The response of the H.E.S.S. array to simulated Cherenkov light is simulated with the detector 
simulation Sim Hessarry which was developed by Bemlohr [2002]. It calculates the PMT re- 
sponse to Cherenkov light for each of the telescopes. Sim Hessarry is a very accurate detector 
simulation which takes into account: 

• the reflector geometry, mirror reflectivity and pointing of the telescopes 

• shadowing by the camera support structure and point spread function, 

• the transmission of the Winston cones in front of the PMTs, 

• the quantum efficiency of the PMTs, 

• the electronic response of the PMTs and 

• the telescope multiplicity requirement of the central trigger. 

Sim Hessarry also simulates the camera's response in different observation modes. The 
standard modes are the on/off and the wobble mode. In the on/off mode, the camera is pointed 
directly towards the source for an on-run and in a source-free direction for an off-run. The 
off-run is used to determine the background. In wobble mode, the camera direction is offset 
from the source direction by the wobble offset is chosen such that the source is still 

enclosed within the FOV. The advantage of this observation mode is that it provides regions 
that are needed for background estimation within a single run. This has the advantage that 
systematic errors on a run by run basis are reduced. For H.E.S.S. observations, the standard 
wobble offset is = ±0.5° in Dec and the corresponding value 6w = ±0.5°/ cos(ZDec) in 
RA. The alternating signs of the wobble offset provide compensation for linear gradients in 
acceptance and hence they reduce systematic errors. 
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5.1.3 Monte Carlo Data 

The Monte Carlo data used in this work is based on a simulation from August 2005. The data 
has been calculated with the desert atmospheric model and partially also with the maritime 
model. These simulations consist of simulated 7-ray showers from point sources. The point 
sources were simulated at zenith angles of 30, 40, 45, 50, 55, 60, 63, 65 and 67° and at 
azimuth angles of and 180° for sources in a northern direction, like the Crab Nebula, and in a 
southern direction, like MSH 15—52. Wobble offsets {0^) were simulated for 0, 0.5, 1, 1.5, 2 
and 2.5°. CT3 was simulated with a reduced efficiency of 8%, as measured in early 2003. At 
each zenith and offset angle, about 3 x 10^ showers reached the sensitive detector area with a 
radius of 1000 m from the center of the H.E.S.S. array. About 10% of these events passed the 
cuts. The energy spectrum of the simulated y-ray showers reaches from 20GeV to lOOTeV, 
with a photon index F = 2. Since the shower simulation consumes most of the computation 
time required, the efficiency of the simulations can be increased if the same simulated shower is 
used multiple times in the detector simulation. The only parameter that is varied is the impact 
position with respect to the telescope. This method can reduce the computational effort by 
orders of magnitudes. The Monte Carlo data is converted and stored in the ROOT file format 
for an easy integration into data analysis. 

Simulations of extended sources have been obtained from the point source Monte Carlo 
simulations by scattering events according to a Gaussian distribution. The width and length of 
the distribution were chosen according to the standard deviation of the extension to be simu- 
lated. Extended simulations were used for the determination of the collection areas. Although 
this method is marginally less accurate than a full simulation it is preferable due to its higher 
computational efficiency. 

5.2 Shower Reconstruction 

The central part of the data analysis is the shower reconstruction. This consists of several steps, 
some of which also require the Monte Carlo data. This section describes the individual steps 
that are taken in order to extract shower information from raw data. 

5.2.1 Camera Calibration 

Calibration is the first step in shower reconstruction, which involves estimating the number of 
photo electrons (p.e.) that have produced a camera event image. A raw image of a camera event 
is shown in Fig. 5.2 (upper left). Pixels below a certain threshold were later removed from 
the data and therefore do not show up in this image. The number of p.e. (A) that hit a PMT is 
calculated from the ADC counts of the PMT's high-gain (A^*^) and low-gain (A^^) channels. If 
A is less than 150 p.e., A^*^ is used. If A is greater than 200 p.e., A}-'^ is used. For intermediated 
values, A is determined from the weighted average of A^^ and A^^ as 

A= (1-e) xADC"^ + exADC^°, (5.1) 

where e = (ADC"*^ - 150)/(200 - 150). 

^^The zenith angle (0) is the angle between the pointing direction of the array and the zenith, i.e. = 90° — 
/Alt, where /Alt is the altitude angle. 
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The amplitudes A^"^ and A^^ are given by 



A''° = ^^^XFF (5.2) 

/e 



for the high-gain channel and by 



ArjpLG _ „LG 

^"^^ = ADC (HG/LG) X FF (5.3) 



e 



for the low-gain channel, where ADC^*^ and ADC'^^ are the recorded ADC counts for the 
high- and low-gain channels, p^^ and p^^ are the pedestals of these channels measured in 

ADC counts. ']^^'" is the conversion factor between ADC counts and p.e. HG/LG is the ampli- 
fication ratio of the high-gain to the low-gain channel and FF is the flat-filed coefficient. These 
parameters are determined with special calibration runs as described in detail by Aharonian 
et al. [2004a] and RoUand [2003]. Since the parameters can vary with time, the calibration has 
to be repeated regularly to guarantee a reliable conversion. 



The Pedestals p^^ and p^*^ 

The pedestal is defined as the mean ADC value recorded in the absence of Cherenkov light. It 
is primarily determined by the night-sky background (NSB) and secondarily by electronic noise 
which has a strong temperature dependence. The pedestals are determined at intervals of 20 s 
from the events of an observation run, where pixels affected by Cherenkov light are excluded. 
Typically only 20 pixels of an image contain Cherenkov light. Fig. 5.2 (top right) shows the 
pedestals during run 20301. Since the width of the pedestal distribution changes with the NSB, 
the pedestal distribution can also be used to determine the level of NSB and to reject noisy 
pixels. 

In addition, every two nights the pedestals are monitored in special electronic pedestal runs, 
which are performed with a closed camera lid. Deviations from the nominal value can clearly 
be measured by these runs. 



The Conversion Factor y^"^ 

The offset of a single p.e. peak from the pedestal determines the conversion factor iy^^^). To 
obtain a histogram of ADC counts that can show the pedestal and the single p.e. peak (Fig. 5.1), 
the camera is illuminated by a faint LED pulser with an intensity of about 1 p.e. per pixel and 
a frequency of 70 Hz. The LED pulser is installed in the camera shelter, where the camera is 
parked during a single p.e. run, in order to be shielded from the NSB while the lid is opened 
to observe the LED pulses. The single p.e. runs have a duration of two minutes and are taken 
every two nights. Since the pedestal ADC counts are Gaussian distributed with a standard 
deviation Op, the ADC counts of a signal of np.e. {n E N) are also Gaussian distributed with a 
standard deviation \fn<3y^ and a mean position P -\- ny^^ . The number of p.e. follows a Poisson 
distribution. Therefore, the conversion factor can be found by the fit to the function / to 
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X /dof 1.060 

Norm 50314 
P -27.31 ±0.26 
16.10 ±0.19 
80.78 ±0.80 
0.46 ±0.01 
1.11 ±0.01 




Figure 5.1: ADC counts of a single 
p.e. calibration run. The first peak 
corresponds to the pedestal counts 
and the second to the single p.e. 
counts. The distance between both 
peaks determines the conversion fac- 
tor 7^'^*^. The height of the single p.e. 
peak is determined by the intensity of 
the LED pulser and obeys the Poisson 
statistic. (Figure taken from Aharo- 
nianetal. [2004a].) 
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the ADC histogram of the single p.e. runs. / is give as 
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,ADC^ 

e y 



nOy, 



(5.4) 
(5.5) 



where x is the number of ADC counts and and K are normalization constants. 



The Amplification Ratio HG/LG 

The amplification ratio (HG/LG) can be determined from a comparison of the ADC counts of 
the high-gain (Ch) and low-gain channels (Cl). Since the pedestals for the high-gain (Pn) and 
low-gain channels (PjJ are also available as described above, 

,5.6) 

LG Cl-Pl 

The amplification ratio is calculated from regular observation runs for all usable pixels with an 
intensity between 15-200 p.e. 



The Flat-Filed Coefficient FF 

Since the calibration of the PMTs is not absolute, differences in the efficiency of the photocath- 
odes or the Winston cones can result in different PMT efficiencies. Deviations with an RMS of 
about 10% have been observed. The FFs compensate for this difference and provide a uniform 
camera response. FFs are determined in special flat-field runs roughly every two days. In a 
flat-field run, the camera is illuminated by an LED flasher which is mounted on the mirror dish. 
At a distance of 15 m from the camera, the flasher provides a homogeneous illumination from 
10-200 p.e. within a solid angle of 10° at the PMTs' peak efficiency in the range of 390 to 
420 nm. The FF of each pixel is determined as the inverse ratio of the pixel amplitude to the 
mean amplitude of all pixels averaged over a run. By definition the mean of the FFs is equal to 
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1. The distribution of the FFs gives an estimate of the uniformity of the camera. The relative 
accuracy after a correction with FFs is < 1%. An image after calibration is shown in Fig. 5.2 
(middle left). 

Identification of Unusable Channels and Broken Pixels 

To obtain a correct calibration, it is also important to eliminate channels that do not provide cor- 
rect information and therefore falsify the image. There are several reasons for unusable chan- 
nels: missing calibration coefficients, synchronization problems with the analog ring sampler's 
memory, high voltage variation in the PMTs due to hardware failures or merely the presence 
of bright stars. If the high-gain and the low-gain channels of a pixel are unusable, the pixel 
is considered "broken". The amount of unusable channels can reach up to 4% per run. The 
number of broken pixels is less. Fig. 5.2 (bottom) shows problematic pixels in a run. 

5.2.2 Muon Calibration 

Camera calibration guarantees images of high quality, but it does not provide an absolute cali- 
bration of intensity since important system components such as the optical components are not 
included. An absolute calibration can be achieved by an analysis of muon rings. This relies 
on the interdependence of the opening angle of the Cherenkov light emitted by muons in the 
atmosphere and muon energy. Since the mirror dish focuses parallel rays onto the same point in 
the focal plane, muons appear as rings in the camera image. The radius of the rings reflects the 
muon's energy and hence the number of emitted Cherenkov photons and the amplitude of the 
muon ring. Muon rings also provide an alternative method for the determination of the flat-field 
coefficients, since neighboring pixels in a muon ring are expected to have a similar amplitude. 
This is explained in Aharonian et al. [2004a]. Details about muon calibration can be found in 
Leroy et al. [2003] and Bolz [2004a]. Muon runs are taken regularly to monitor the H.E.S.S.' 
photon efficiency, which decreases by about 5% percent per year (Bolz [2004a, pg. 89]). The 
degradation is partially compensated for by regular hardware maintenance and upgrades every 
few months, keeping the deviations from the nominal values sufficiently low. Nevertheless, cor- 
rections for this absolute efficiency have been developed by Bruno Khelifi and Conor Masterson 
[2005] and can be implemented in the future. 

5.2.3 Image Cleaning 

The purpose of image cleaning is to remove the pixels from an image which do not represent 
Cherenkov light from the shower. The cleaning is realized with a tail cut algorithm that de- 
termines which pixels will be kept and which will be removed from the image. The tail cut 
algorithm has two cut parameters (l,h) — the low and the high thresholds of the pixel ampli- 
tude. The standard values are /=5p.e. and /?=10p.e. The tail cuts remove all pixels with an 
amplitude below / and keep all pixels with an amplitude above h, if they have a neighboring 
pixel of at least /. They will also keep any pixels within the range of / and h if they have a 
neighboring pixel with an amplitude exceeding h. The result is a cleaned shower image, as 
shown in Fig. 5.2 (middle right). 
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Figure 5.2: Camera images of CT4, run 
20301, event 4923 at different steps in the 
analysis: the zero suppressed high-gain ADC 
counts (top left), the pedestal ADC counts 
(top right), the corresponding pixel intensi- 
ties in p.e. after calibration (middle left), the 
cleaned image after (5,10) tail cuts (middle 
right) and the unusable channels and broken 
pixels during the run (bottom). The Hillas el- 
hpse with its center of gravity is indicated in 
each image. 
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5.2.4 Image Parameterization according to Hillas 

The standard analysis uses Hillas parameters, which are derived from cleaned camera images. 
Hillas parameters were first introduced by Hillas [1985]. They are the moments of an image 
as described in App. B. Their geometric representation is shown in Fig. 5.3. They provide a 
powerful means for shower reconstruction and the rejection of background events. The relevant 
parameters for analysis are the image's amplitude (M, Eqn. B.2), local distance (LD, Eqn. B.5), 
width {W, Eqn. B.17), length (L, Eqn. B.17) and angle Eqn. B.20). 



5.2.5 Geometric Reconstruction 

Stereoscopic lACTs allow geometric shower reconstruction (Aharonian et al. [1997], Hofmann 
et al. [1999]) because the shower direction and impact point are uniquely defined by the shower 
images of two telescopes. With more than two telescope images, shower geometry is even 
over-defined and reconstruction accuracy increases significantly. The increase results from the 
fact that more than two telescopes build an array, which guarantees that the majority of shower 
images will be taken at a favorable angle <^ 1 80° and distance < 50 m. The actual reconstruction 
is based on the geometric optics of lACTs. Reconstruction of the shower direction and the core 
position, i.e. the shower impact point at the ground, is done in two different coordinate systems 
that are associated with the telescope array (cf. Gillessen [2004], Ergin [2005]). 

Reconstruction of the Shower Direction 

Shower direction is reconstructed in the nominal system. The nominal system is a two-dimen- 
sional coordinate system perpendicular to the line of sight of the array. It represents the camera 
coordinate system in an imaginary nominal focal plane with a focal length of 1 m. Distances 
in this plane are measured in radians and correspond directly to sky coordinates. Moreover, in 
the nominal system the projected shower direction is found along the elongated major axis of 
the Hillas ellipse. Hence the shower direction is determined by the intersection of the elongated 
major axis of a camera images of an event. In a perfect reconstruction the lines intersect in one 




Figure 5.3: A geometric representation of 
Hillas parameters. The ellipse represents 
a y-ray shower image. (Sketch taken from 
Ergin [2005, pg. 80].) 
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point {x,y). However, due to limited reconstruction accuracy, the pairwise intersection points 
{x,y)ij of the lines // and Ij = 1,2,3,4) differ as shown in Fig. 5.4. So a unique definition 
of the shower direction {x,y) is required. Different definitions are possible. The H.E.S.S. stan- 
dard analysis has adopted the following definition which is based on the weighted sum of each 
intersection. The intersection point is 



with the weight 



W; 



sin(0/-0/) 



(m; + iaMtjw: + ■ 



(5.7) 



(5.8) 



where Wij is an empirical weight. The weight is defined by three terms accounting for three 
different quality criteria. The term in the numerator favors the intersection of orthogonal lines, 
which have a smaller error factor than parallel lines. The I A term in the denominator puts more 
weight on images of higher image amplitude which provide a higher accuracy. The L/W term 
suppresses images with a small length to width ratio, where the major axis is less well-defined. 
Determining the shower direction is necessary for the generation of sky maps and subsequent 
statistical analysis. 



Figure 5.4: Reconstruction of shower direc- 
tion from four shower images in the nomi- 
nal system. The intersection of the elongated 
major axis determines the shower direction. 
(Figure taken from Ergin [2005, pg. 85].) 




Reconstruction of the Shower Core Position 

The shower core position is determined in a tilted system, similar to the shower direction. The 
tilted system is a three-dimensional coordinate system perpendicular to the line of sight of the 
array. Each camera image has a distinct position determined by the center of the telescope. 
Since the shower direction is found along the major axis of the Hillas ellipse, the intersection 
of the elongated major axis of all images in the tilted system defines the shower core position 
(Fig. 5.5). Transformation to the ground system provides the shower impact point at the ground. 
The core position is also defined by the weighted sum of the pairwise intersection {Xjy)ij given 
by Eqn. 5.7, but with different empirical weights. Either 

= |sin(<^,--<^^)| (5.9) 
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or 



Wij 



(5.10) 



is used. The difference is minor. For the analyses presented here, Eqn. 5.9 has been used. The 
shower core position is needed to calculate the impact distance, which is the distance of the 
core position to the telescope. This parameter is important for the energy reconstruction of the 
shower. 




Figure 5.5: Reconstruction of the shower 
core position from four camera images in 
the tilted system. The intersection of the 
elongated major axis determines the shower 
core position. The camera images are greatly 
magnified in comparison to their distance. 
(Figure taken from Ergin [2005, pg. 86].) 



5.2.6 Angular Resolution and Point Spread Function (PSF) 

The angular resolution of H.E.S.S. is determined by the point spread function (PSF), which 
is given by the accuracy of the geometric reconstruction of the shower direction. The PSF 
can be determined from Monte Carlo simulations. The H.E.S.S. PSF has been simulated and 
parameterized for the relevant configurations. The parameterization is used where information 
about the PSF is required, e.g. to indicated the PSF in sky maps. Fig. 5.6 shows sample 
distributions of the reconstructed shower direction versus the square of the angular distance 
(0) from the true value. These 0^-distributions are obtained from point source Monte Carlo 
simulations for different zenith angles. The distributions in Fig. 5.6 have been normalized 
such that the integral is one. The parameterization (PSF) is given by the sum of two Gaussian 
functions as 

-e2\ . /-02- 



PSF{G] 



(5.11) 



with the standard deviations Oi and 02, the relative amplitude A and the normalization con- 
stant A. The fit functions are shown by the solid lines of the same color. Small deviations can 
be seen near the origin. The representation versus 6^ should not be confused with a profile of 
the PSF which is linear in 6. The width of the PSF can be expressed by a single value: the 68% 
containment radius (r68%). It is defined as the radius of a circle at the signal regions, which 
contains 68% of events in case of a point source. 
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The PSF depends on the zenith angle (0), wobble offset azimuthal orientation (Az) 
during observations (due to the magnetic field of the earth) and image amplitude cut (M). For 
example, the change of the image amplitude cut from 80 to 400 p.e. results in a reduction of 
the width of the PSF by ~50%. The PSF also depends on the y-ray energy spectrum, but this 
dependence is small and negligible for the data discussed in this work. A summary of different 
PSF fit parameters and r(,g% is given in Tbl. 5.1. 



HESS PSF for Different Zenith Angles 
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Figure 5.6: 0^ -distributions of the H.E.S.S. 

PSF as determined from point source simu- 
lations at different zenith angles (0), d„ = 
0.5°, south orientation and lA > 80p.e. The 
distributions are fit with the sum of two 
Gaussian functions (solid lines). The con- 
tainment radii of 68% are indicated by the 
dashed vertical lines. 



Table 5.1: Fit parameters and 68% containment radius (rgg^) of the H.E.S.S. PSF model for different 
zenith angles (0), wobble offsets (6^), azimuthal orientations (Az) and image amplitude cuts (M). 
The relative error of the parameters Oi, G2, Ai-ei and r68% are 0.5%, 0.5%, 2% and 1% respectively. 
The parameterizations for different configurations have been numbered to simplify a later reference. 



No. ®[°],e4°],Az,IA[p.e.] Are! 



02 



1 30, 0.5, south, 80 0.0601 0.0477 0.111 0.112 

2 40, 0.0, south, 80 0.175 0.0481 0.117 0.121 

3 40, 0.5, south, 80 0.170 0.0477 0.117 0.120 

4 40, 1.0, south, 80 0.166 0.0472 0.116 0.118 

5 40, 0.5, south, 400 0.0655 0.0333 0.0785 0.0632 

6 50, 0.5, south, 80 0.206 0.0489 0.125 0.139 

7 50, 0.5, north, 80 0.174 0.0458 0.121 0.129 



5.2.7 Energy Reconstruction 

Energy reconstruction relies on the dependence of image amplitude (M) on 7-ray energy (E), 
impact distance (b) and zenith angle (0) . This dependence of IA{E, b, 0) can be used to calcu- 
late the energy E{IA,b, 0) as a function of image amplitude, impact distance and zenith angle. 
This is accomplished with Monte Carlo simulations. For each telescope and certain zenith an- 
gles, lookup tables of the energy as a function of I A and true b are created. For a shower recorded 
by a telescope with the image amplitude lAj^i, the impact parameter bj^i and the zenith angle 0, 
the energy ^xeK^^Teb^Tei,®) can be obtained from the corresponding lookup table. Values for 
intermediate zenith angles are found by linear interpolation in cos(0). The shower energy is de- 
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Energy Lookup Table (0=40°, CT1 ) 
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Figure 5.7: Energy lookup table of CT 1 for the 
zenith angle of 40°. The mean reconstructed en- 
ergy is represented as a function of the image 
amplitude and the simulated impact distance. 




Figure 5.8: Distribution of the relative error 
in the reconstructed energy per event for simu- 
lated y-rays. (Figure taken from Aharonian et al. 
[2006a].) 



termined as the arithmetic mean of the energies which have been determined for each telescope. 
Fig. 5.7 shows the lookup table of CT 1 for a zenith angle of 40°. 

The relative error (A^) for a simulated y-ray shower with Monte Carlo true energy Etme and 
reconstructed energy ^reco is defined as Ae = (^reco — £^true)/£^true- Above the energy threshold, 
A^ has a standard deviation (cr^) of ~15% equivalent to the energy resolution of a single event. 
Fig. 5.8 shows Oe as a function of ^tme for simulations with a power law energy spectrum 
(r = 2.6) above 440 GeV at a zenith angle of 50°. The RMS is 16% and the width of the 
Gaussian fit 14%. 



5.3 Statistical Methods 



Statistical methods are necessary to find and characterize the y-ray signals in the data. The 
relevant techniques are discussed here and are demonstrated with data from the Crab Nebula. 



5.3.1 Run Selection and Quality Criteria 

Before starting a data analysis it is important to verify the data quality. For example, bad 
weather or technical problems are a few reasons why data of low quality could be a part of a 
data set. Since this data mainly adds background to the analysis and produces systematic errors, 
this data should be excluded from the analysis. The following empirical quality criteria have 
been developed in order to decide whether an observation run should be included in an analysis 
or not. 



Trigger Rates 

The system trigger rate (R) is a very sensitive measure for the quality of atmospheric conditions. 
Since the hadronic background is very constant, deviations from the typical trigger rate can 
indicate variations of the transmissibility of the atmosphere (Fig. 5.9). Since such data cannot 
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be handled in the standard analysis, these runs have to be excluded. The individual trigger 
criteria for including a run are: 

• Root mean square (RMS) ofR< 30% 

• A relative change ofR< 30% 

• RMS of R after subtraction of the fit of by a straight line < 10% 

• R> ^thres- The threshold Rthiesi®j) is derived from reference observations taking into 
account the zenith angle (Funk et al. [2004]) and exponential efficiency loss over time (t) 
(H.E.S.S. collaboration [2001]). 
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Figure 5.9: Trigger rate of run 20418. The sta- 
bility is typical for a run included in an analysis. 
The slight slope reflects the zenith angle depen- 
dence of the trigger rate. 



Broken Pixels 

Broken pixels and unusable channels (Sec. 5.2.1) can lead to systematic errors in reconstruction. 
The criteria for including a run consist of: 

• Number of broken pixels during a run for any camera < 120 

• Number of pixels with only the high voltage turned off < 50. 



Dropped Events 

Sometimes not all of the triggered events can be recorded due to technical problems. If the 
number of dropped events exceeds 5%, the problem is considered as severe and the run is 
excluded from the analysis. 



Tracking Accuracy 

To maintain the system's pointing accuracy of 20 arc seconds, runs with a standard deviation of 
tracking more than 10 arc seconds are excluded from the analysis. 
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5.3.2 Event Selection 

The selection of events with certain characteristics is useful to separate y-ray showers from the 
background. The event selection is accomplished with special selection cuts. These cuts reduce 
only a small fraction of the 7-ray showers, but a large fraction of the background. An important 
aspect of the event selection is 7 / hadron separation. For example, while the typical system 
trigger rate is ~300 Hz, the actual y-ray rate is less than ~0. 1 Hz. However, with the right 
selection cuts this small signal to noise ratio can be increased by a few orders of magnitude. A 
short description of useful selection cuts is given below. 

Quality Cuts 

The main purpose of quality cuts is to provide a reliable shower reconstruction and an efficient 
background reduction. Three such cuts are used: the cuts on the trigger multiplicity, those on 
the image amplitude and those on the local distance. 

The trigger multiplicity cut selects all events with a minimum number of triggered tele- 
scopes. Its default value is two providing stereoscopic events, only. The image amplitude cut 
improves the accuracy of shower reconstruction. It removes images of small image amplitude 
that do not have a well-defined major axis and therefore are difficult to reconstruct. The standard 
value is 80p.e. The local distance (LD) cut controls the distance between the center of gravity 
of an image and the camera center in the nominal system. A standard value of 2° guarantees 
that the majority of the images are not truncated by the camera edges, which would falsify the 
reconstruction. The values of the quality cuts are listed in Tbl. 5.2 and 5.3. 

Mean Reduced Scaled Cuts 

The Hillas parameters of width (W) and length (L) provide efficient means for 7 / hadron 
separation and background suppression. The separation is achieved through the different shapes 
of 7 and hadronic shower images (cf. Sec. 3.1). However, a separation cannot be achieved with 
a cut on the width or length alone, since width and length also depend on the image amplitude 
(I A), the impact distance (b) and the zenith angle (0). Nevertheless, Monte Carlo simulations 
can reproduce this complex dependence and provide this information in lookup tables with the 
mean width (< W >) and length (< L >) and the corresponding standard deviation (ow and Ol) 
as a function of the natural logarithm of the image amplitude and the impact parameter. These 
lookup tables are calculated for different zenith angles where intermediate values are obtained 
by linear interpolation in cos(0) . Fig. 5.10 and 5.11 show the lookup tables for a zenith angle of 
40° with a wobble offset of 0.5°. With these lookup tables, one can calculate the mean reduced 
scaled width (MRSW) which is defined as 



where A'xei is the total number of telescopes. The mean reduced scaled length (MRSL) is de- 
fined analogously. Hence the MRSW and MRSL express the deviation of width and length of 
a shower image from its expectation value. The deviation is measured in units of standard 
deviations. Events with small MRSW and MRSL parameters indicate high similarity to 7-ray 
air showers and vice versa. Therefore, cuts on high MRSW and MRSL parameters reduce the 
number of background events. Fig. 5.12 and 5.13 show the MRSW and MRSL distributions 



MRSW= 




(5.12) 
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Mean Width Lookup Table (0=40°) 




ln(IA/p.e.) 



Figure 5.10: Lookup table of mean width for a 
zenith angle of 40°. 



Ow Lookup Table (0=40°) 




ln(IA/p.e.) 



Figure 5.11: Lookup table of the standard devi- 
ation of the width for a zenith angle of 40°. 



of simulated y-ray showers and background showers from source free observations at a zenith 
angle of ~ 40° and at a wobble offset of 0.5°. The y-ray showers are approximately Gaussian- 
distributed around zero with a standard deviation of one. The background events are asymmet- 
rically distributed, with a tail extending to high values. The mean reduced scaled standard cuts 
are indicated by the vertical lines. They remove a large fraction of background events while 
keeping most y-events, resulting in an effective 7/ hadron separation. The mean reduced scaled 
cuts for different cut configurations are listed in Tbl. 5.2 and 5.3. 



MRSW (0=40°, 9,„=0.5°) 




MRSW 

Figure 5.12: Mean reduced scaled width dis- 
tribution for background and Monte Carlo 7-ray 
data at a zenith angle of 40° and a wobble offset 
of 0.5°. 



MRSL (0=40°, 9w=0.5°) 




MRSL 

Figure 5.13: Mean reduced scaled length dis- 
tribution for background and Monte Carlo 7-ray 
data at a zenith angle of 40° and a wobble offset 
of 0.5°. 



5.3.3 Background Models 

Although a large fraction of background events can be reduced with cuts, some 7-ray-like back- 
ground events cannot be distinguished and separated from 7-ray events and contaminate the 
true 7-ray counts. Nevertheless, the true strength of a 7-ray signal can still be determined if 
this contribution from the background is taken into account. The strength of the background is 
determined from background regions, which are defined by a background model. The signal 
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and background regions are referred to as ON- and OFF-regions. Various background models 
exist. The region- and the ring-background models are used in the standard analysis. Both 
use a circular ON-region with a radius 6 and a size proportional to 6^. The OFF-regions are 
located at some distance from the ON-region where virtually no y-ray showers are observed. 
The H.E.S.S. FOV of 5° is large enough to allow for simultaneous ON- and OFF-regions in the 
same run. This is an advantage, since systematic errors on a run by run basis are reduced. Also, 
bright stars can contribute to systematic errors. As shown by Puehlhofer [2004], stars brighter 
than a magnitude of about six can cause dips in acceptance at the star position. The dips are 
presumably caused by reduced image amplitudes from switched-off pixels. To reduce these 
systematic errors, stars in the ON- and OFF-regions should be avoided. The advantage of the 
region-background model is its compensation for the radial gradient of acceptance in the FOV, 
which makes it the first choice for spectral analysis and robust against systematic effects. On 
the other hand, the advantage of the ring-background model is its provision of sky maps. 

The Region-Background Model 

The region-background model can be applied if the data was taken in wobble mode. The ON- 
region is chosen at the source region. The distance to the center of the FOV is given by the 
wobble offset (d^). It radius (9) is limited by this distance {0 < 0^)- The OFF-regions are 
placed on a circle given by the center of the FOV and a radius equal to the wobble offset. They 
are preferably arranged symmetrically to the ON-region with respect to the center of the FOV. 
This choice provides a very similar acceptance for the ON- and OFF-regions independent of the 
radial gradient of the acceptance. The number of OFF-regions may vary. More OFF-regions 
provide higher background statistics. The maximum number of OFF-regions is limited by their 
radius and wobble offset. In addition, a certain distance to the ON-region has to be preserved to 
avoid a contamination of the OFF-regions with 7-ray events. Fig. 5.14 illustrates the geometric 
layout for two runs with a wobble offset of ±0.5° in Dec and seven OFF-regions each. The 
distance to the ON-region is given by the 95% containment radius of the source. The number 
of OFF-regions (A'oFF-regions) determines the normalization constant a. 



1 



(5.13) 



a = 



NoFF— regions 



FOV wobble -0.5° 



FOV wobble +0.5° 




Figure 5.14: Illustration of the 
region-background model applied to 
two runs with wobble offsets of 
±0.5° in Dec and seven OFF-regions 
(grey) each. The ON-region (black) 
is located in the center. The contain- 
ment radius of 95% defines an area 
where no OFF-region is contained. 
(Figure taken from Piihlhofer [2001, 
pg. 149].) 



regions 
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-1.0 
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The Ring-Background Model 

The ring-background model can be applied to data regardless of its observation mode. The 
ON-region is circular and has a radius (6). The OFF-region is a ring which encloses the ON- 
region. The area of the ring is about seven times the area of the ON-region. The radius of the 
ring is chosen close to the signal region but not so close as to be significantly contaminated 
by the signal in the ON-region. A typical value for point sources is 0.5°. To obtain sky maps, 
this background model is applied to each bin in a sky map. This is illustrated in Fig. 5.15 for 
six different sky positions. The results obtained at the source position provide the statistics of 
the source. Sectors which contain y-ray regions are excluded from the ring-background. For 
example, the ON-region is always excluded. However, it is not only the ratio between the ON- 
and OFF-region that determines the normalization constant a, but it is also the acceptance. 
Therefore a is determined as 

" = 7 e(xy)dxdy ' withx,y^R,^ (5.14) 

JOFF-region I-*' Jjaxay 

where x and y are the sky coordinates and e{x,y) is the acceptance for the data set. £{x,y) is 
given as 

e{x,y) = — £rm{x,y)tiun, (5.15) 

i^runs ^run nins 

where (?run) is the live-time and £mn{x-,y) is the acceptances for each run. The live-time is the 
dead-time-corrected observation time and emn {x, y) is obtained from a fit to the radial acceptance 
profile. 

Due to the integration over the ON-region, the resolution of a sky map provided by the 
ring-background is limited by the size of the ON-region. The resolution can be increased if the 
ON-region is reduced to the size of a bin. The excess maps discussed in this work are of this 
type. However, alpha is still calculated with the original ON-region. 

A disadvantage of the ring-background model is the energy dependence of the acceptance, 
which makes spectral analysis more difficult. Also, for technical reasons the calculations are 
obtained from count maps where some accuracy is lost through the binning. This is the rea- 
son why small differences from the region-background model can be observed. However, the 
difference is insignificant and negligible. 



FOV wobble -0.52 




FOV wobble +0.5° 



• ^background 
■' control 
ring 

^exclusion 
region 



+ 



+ 



+ 



-1.0 -0.5 +0.5 +1.0 
declination - source declination [deg] 



Figure 5.15: Illustration of the ring- 
background model. At each sky posi- 
tion a pair of signal- (black) and ring- 
regions (grey) is defined. Here only 
six pairs of rings are shown. Regions 
containing /radiation, e.g. the signal 
region, are excluded from the ring- 
background. 
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5.3.4 Excess and Significance 

With a background model one can calculate signal statistics. The number of true y-ray events is 
given by the y-ray excess (Nj), which is determined as 

Nyr = A^oN - aA^OFF, (5.16) 

where A^qn and A^off are the number of events in the ON- and OFF-region and a is the nor- 
malization constant given by the background model. The error of excess {ANj) is obtained 
according to Gaussian error propagation as 



ANy - -^/AA^^N + («AA^off)2, (5.17) 

where AA^qn = V^ON and AA^qit = V^off provide good approximations for statistical errors. 
The y-ray rate is obtained by division with the live-time (t), i.e. 

rate = -'^. (5.18) 
t 

The significance of the y-ray signal in the ON-region according to Li and Ma [1983] is 



S=V2 



AToKln ( ) + AToppln ( (1 + a) ^"^^ ] 
V « A^on+A^off; V ^No^ + NoffJ 



1/2 



(5.19) 



In y-ray astronomy, the significance is an important value since until recent years only a few 
sources had been known, but a lot are now being detected. The detection of a source is only 
claimed if it is observed with a significance of at least five standard deviations (S > 5). A 5a 
detection implies the probability for the signal being caused by a statistical fluctuation of less 
than ~ 6 X 10"^. 

The significance depends on the cuts used in the analysis, which can be optimized to give the 
maximum significance. The optimized cut configurations of the standard analysis are presented 
in Aharonian et al. [2006a] and are listed in Tbl. 5.3. They vary with the strength and the 
spectrum of the source. The extended cut configuration is optimized for a flux of 10% of the 
Crab Nebula with a similar spectrum from an extended source. The set of hard cuts is optimized 
for sources with a flux of 1% of the Crab flux and a photon index (F) of 2.0. The loose cuts are 
optimized for strong sources similar to the Crab Nebula with a F of 3.0. The cuts used here are 
based on these cut configurations and differ only in the 0^-cut, which has been adapted for an 
extended source ((7^ = 0.1°). 



Table 5.2: Cuts parameters which are used in all cut configurations of the standard analysis. 



Trigger Mult. 


LD 


MRSL 


MRSL 


MRSW 


[No. oftels] 


[deg.] 


Min. 


Max. 


Min. 


2 


2.0 


-2.0 


2.0 


-2.0 
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Table 5.3: Different cut configurations of the standard analysis. 6^^^ denotes the standard cut values 
and 0^ the modified values which have been adjusted for the analysis of MSH 15—52. Additional 
cuts parameters which are common to all configurations are listed in Tbl. 5.2. 



Cut-Configuration 


MRSW min. 


I A min. 




02 






[p.e.] 


[deg.2] 


[deg.2] 


extended 


0.9 


80 


0.16 


0.09 


hard 


0.7 


200 


0.01 


0.09 


loose 


1.2 


40 


0.04 


0.09 



5.3.5 Source Position and Size 

The y-ray excess map obtained from the ring-background model can be fit with the two-dimen- 
sional Gaussian function G(Tr,cry,a(),5(),<B,A'(^' ^) determine the position, size and orientation of 
a 7-ray signal. G is obtained through the convolution of a two-dimensional Gaussian with the 
PSF as described in App. C. This way it is possible to distinguish between the influence of the 
PSF and the true i.e. intrinsic width of the 7-ray source. The parameterization of the H.E.S.S. 
PSF according to Eqn. 5.11 is determined with Monte Carlo simulations (cf. Sec. 5.2.6) and 
summarized in Tbl. 5.1. The parameters Uq and 5o represent the fit position in RA and Dec. 
The parameters Ox and Oy are the standard deviations of the intrinsic width and length of the 
y-ray signal. (O is the angle between the major axis of the fit function and the RA-axis. 

5.3.6 Analysis of Data from the Crab Nebula 

To demonstrate and verify the methods described above, they are applied to H.E.S.S. data from 
the Crab Nebula taken from January 25, 2004 to March 4, 2005. The data was taken in wobble 
mode with offsets of 0vv = ±0.5° in RA and Dec at (J2000) (15i^l4™27% -59°16'18"). The 
mean zenith angle of the observations is 50.0°. The complete run list of the 30 runs with a 
live-time of 12.4 h after run selection is given in Tbl. E.l of App. E. The analysis was carried 
out with extended cuts (Tbl. 5.3). 

Fig. 5.16 and Fig. 5.17 show the excess and the significance map as obtained with the ring- 
background model. The configuration of the ring-background and the region-background model 
are indicated. The ON- and OFF-regions have been chosen to not contain stars brighter than 
a magnitude of six. The ON-region has a radius of 0.3° explaining the broad peak of similar 
radius in the significance map. The inner and outer radii of the ring-background are 0.725° and 
1 .075° respectively. Also the region-background model was applied using one OFF-region for 
each of the four wobble offsets. The statistics of both models are summarized in Tbl. 5.4. A 
very significant y-ray signal is detected at the position of the Crab Nebula. 

The results from a fit of the excess map with the Gaussian fit function G(7^,(7y,ao,^,fi),Af ^) 
are shown in Tbl. 5.5. The fit position is found close to the pulsar position as shown in the excess 
map of Fig. 5.18. The intrinsic width and length are close to zero. The data and the fit function 
along two orthogonal axes are shown in Fig. 5.19. The extension of the excess can be explained 
by the component of the PSF alone. Taking the system's pointing accuracy of 20" into account, 
the best fit position of the Crab Nebula is (J2000) (5'^34"i30.^6±4.^5stat± l?4syst,22°l'9;'8± 
3!'4stat±20',yst)- For comparison, the position of the Crab Pulsar (PSRB053 1-1-21) as deter- 
mined from radio observations is (5*^34^31.^97, 22°0'52f'07). 
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05h40m 05h30m 

RA (hours) 



Figure 5.16: Excess map of the Crab Neb- 
ula showing the configuration of the ring- 
background model as applied in the analysis. 
Stars are indicated with their magnitudes. 




RA (hours) 



Figure 5.17: Significance map of the Crab Neb- 
ula showing the configuration of the region- 
background model as applied in the analysis. 
Stars are indicated with their magnitudes. 



Table 5.4: Signal statistics as obtained from the ring-background and the region-background models 
for 30 runs from the Crab Nebula with a total live-time of 12.4 h. The marginally different event 
statistics in the ring-background model are explained by binning of the data. 



Ring-background Region-background 



A^ON 


11479 


11512 


NoFF 


25485 


4380 


a 


0.180 


1 


Ny 


6899 ± 111 


7132 ± 126 


S[a] 


76.3 


57.6 


signal/noise 


21.7 
1.50 


16.4 
1.63 


rate [min^] 


9.3 ±0.12 


9.62 ±0.17 



Table 5.5: Fit parameters of the two-dimensional Gaussian function Gg- (j,,,ceo,5o fi).A' 

{a, 5) (App. C) 

for the y-ray excess map of the Crab Nebula. The parameterization of the PSF is found in Tbl. 5.1. 



Parameter 


Value 


Position RA {Oq) 
Position Dec (5o) 
Length (c7^) 
Width (Oy) 
PSF parameterization 


5'^34'"30.^6±4.^5, (83.627° ±0.001°) 
22°1'9.84" ± 3'.'4, (22.019° ± 0.001°) 
0.3' ±0.3', (0.005° ±0.005°) 
0'±0.3', (0° ±0.005°) 
no. 7 (Tbl. 5.1) 
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Figure 5.18: Excess map of the Crab Nebula 
with the contour hnes (white) of the Gaussian 
fit function. The fit position of the centroid 
and the systematic errors are represented by the 
black cross. The position of the Crab Pulsar 
PSRB0531+21 is indicated by the black circle. 
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Figure 5.19: Slices along two orthogonal axes 
of the 7-ray excess. The Gaussian fit function, 
the PSF and the intrinsic width and length are 
indicated. Within statistical errors the intrinsic 
width and length are zero (cf. Tbl. 5.5). 



5.4 Spectroscopy 

Spectroscopy is an important aspect in astronomy because it can provide new information about 
an astrophysical object and the processes of a source region. If 4> is the flux from the source, 
then its energy spectrum is given by the differential flux Often the energy spectrum obeys 
a power law of the form 

where ^ixeV^ T and E denote the differential flux at 1 TeV, the photon index and the energy, 
respectively. The integral flux above 1 TeV can be obtained by integration of the differential 
flux as 

4>(£>lTeV)=r ^iTevf-TT^) dE = ^TcY. (5.21) 

JiTev viTevy 1-r 

The integral flux is used for comparison of the flux measured by experiments with different 
energy thresholds. The error A4>(£' > ITeV) is obtained by Gaussian error propagation.^^ If 



According to Gaussian error propagation, the error of the function g{x,y) is given by 



where a stands for the standard error and cov(jic,3') for the covariance between x and y. Here g, x and y have to be 
substituted by 0iTeV and F, respectively. 
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CiTeV; c^r and cov(aiTeV)r) denote the error of the differential flux, the error of the photon 
index and the covariance, then 



o^{E > ITeV) = 




2 I ( '^ITeV 
<l>lTeV^ \ (1 -r)2 



01TeV 



cov((^iTeV,r)> (5.23) 



where the standard errors and the covariance are obtained from a ;^^-fit to the energy spectrum. 

The energy spectrum can be determined from the number of 7-ray excess events (Ny) pro- 
vided by the background model. If the detector's effective area (A) is known, the differential 
flux is given by 



dE 



1 dNy 

WdE 



1 

1^ 



NON 1 A^OFF 1 



i=l 



■ 1 A, 



(5.24) 



Here t is the live-time of the observation, A^qn and A^off are the number of events in the ON- 
and OFF-regions and A; is the effective area for an event /. The normalization constant a as well 
as the ON- and OFF-regions are determined by the background model. Since the usage of the 
ring-background model is difficult in spectral analysis, in this work only the region-background 
model is used for spectroscopy. Again, the error of the flux o (^) is found by Gaussian error 
propagation for Eqn. 5.24 as 



d^ 
'dE 



tAE 



1 

7ae^ 



NoN ^ ^ NoFF ^ ^ 

^ 42 Af ^' f^^ A] A4 ^1 



(5.25) 



1=1 



Here C7(A,) is the error of the effective area for an event /. The covariance is zero, since the 
error of the event statistics and the effective area are uncorrelated. Eqn. 5.25 takes the statistical 
error of the events and the effective area for each event into account. 

In the analysis, a histogram with differential flux and a logarithmic binning is filled with 
the excess given by Eqn. 5.24. The bin size is chosen by an empiric rule, according to the 
significance of the source. Six bins per decade are typical for strong sources. The spectral 
parameters 0iTeV and T are obtained from a least-square fit to the flux. The fit range covers all 
bins from the first bin which is above the safe energy threshold up to the last significant bin with 
a relative error of less than 95%, i.e. two standard deviations. The bin errors represent the 68% 
Feldman-Cousins' confidence interval which is discussed in App D. These asymmetric errors 
are also taken into account in the X^-^^- If the fit value for a bin is greater (or smaller) than the 
bin content, then the positive (or negative) error is chosen for the corresponding bin error of the 
X^ fit. 



5.4.1 Effective Area 

The effective area, also known as collection area, of a detector is the corresponding area of 
an imaginary detector with an efficiency of 100% that would detect the same event rate. The 
effective area (A) of the H.E.S.S. detector with a sensitive area {H) of about 5 x 10^ m^ is a 
function of the energy {E), the zenith angle (0), the azimuth angle (^) due to the magnetic 
field of the earth and the wobble offset (0^)- A(£',0, 0vv) is determined by Monte Carlo 
simulations from the number of detected events past cuts A^(£', 0, 0, 0^) and the total number 
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of simulated events Nioia\{E, 0, 0, 9^) as 



Motal(£,0,<^,0w)' 



(5.26) 



Since A is dependent on the cuts, A differs for each cut configuration. To reduce systematic 
errors from the bias of the energy reconstruction discussed in the next section, A is determined 
as a function of the reconstructed (£reco 

) instead of the true shower energy (^tme) • After A is 
determined, it is fit by the empirical fit function 



f{E) = po ■ + p2E^ + p3E^ + p^E^ + p^E + p(,. 



(5.27) 



^(£'reco) has been calculated for the Monte Carlo files of Sec. 5.1.3. Intermediate values are 
determined by linear interpolation in /, cos(0) and the wobble offset Ow 

Fig. 5.20 shows the effective area of the extended cuts (Aext) for different zenith angles. 
Fig. 5.21 shows Aext in comparison with different effective areas. Atme is the effective area as 
a function of the true instead of the reconstructed energy. The difference is apparent at low 
energies. Agxt is the effective area for the extended cuts and an extended source with a standard 
deviation in width {Ow) and length (C/) of 0.04° and 0.11° respectively. Since its difference 
from Aext is only minor, it is not apparent in the figure. 




E [TeV] 



Figure 5.20: Effective area for different zenith 
angles for the extended cuts used in the analysis 
of Crab Nebula data. 




E[TeV] 



Figure 5.21: Effective area for different cut 
configurations as described in the text. The dif- 
ference between the effective area of true and 
reconstructed energy is apparent. 



5.4.2 Energy Bias and Systematic Errors 

The accuracy of spectral reconstruction is determined by the energy reconstruction. The mean 
error of the reconstructed and the true energy (A^ = (^reco — £^tme)/£^true) is shown in Fig. 5.22 
as determined from simulations. A bias of A^ towards both ends of the H.E.S.S. energy range 
can be seen. The reason for this is that it is a bias in the event selection which favors events 
within the H.E.S.S. energy range and suppresses events at the ends. To avoid the resulting 
systematic errors in A and the energy spectrum, A is calculated as a function of reconstructed 
instead of the true Monte Carlo energy. A disadvantage is that A becomes dependent on the 
energy spectrum of the simulations. If this dependence was strong, each of the sources would 
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Figure 5.22: Bias in the energy reconstruction 
{Ereco — Etme) / Etrue and Safe energy threshold 
for various zenith angles. The bias above the 
safe energy threshold is less than 10%. 
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Figure 5.23: An alternative definition of the en- 
ergy threshold as the peak value of the differ- 
ential event rate dR/dE past cuts. For compar- 
ision, the energy threshold is indicated for the 
same zenith angles as in Fig. 5.22. The safe 
threshold is slightly higher. 



require their own A which could be determined in an iterative approximation of the simulated 
to the true spectrum, as described in HEGRA collaboration, F. A. Aharonian et al. [1999]. For- 
tunately this is not required in most cases. Another problem is the steepening of A near the 
energy threshold, which makes this energy range very sensitive to systematic errors. To avoid 
this problem, the energy range for spectroscopy is restricted to the range with an energy bias of 
less than 10%, which defines the safe energy threshold. Fig. 5.22 shows the safe energy thresh- 
old for different zenith angles. The safe energy threshold is a bit higher than the maximum of 
the differential event rate dR/dE past cuts, which is commonly defined as the energy threshold 
(Hofmann et al. [1999]). Using Eqn. 5.20, the energy threshold is given by 



d^ 



dR 

Ie 



f E \ 



(5.28) 



Fig. 5.23 shows this alternative definition of the energy threshold to the safe energy threshold. 
It provides similar results but it is not used in the standard analysis. 

The systematic error due to differences in the spectrum between data and simulations was 
studied in detail in Aharonian et al. [2006a] and Berge [2006]. It has been verified that spectra 
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Figure 5.24: The ratio of the reconstructed to 
the true effective area per energy bin, for photon 
indices from 1.1 to 3.2, based on Monte Carlo 
simulations at a zenith angle of 45° and with a 
photon index of 2.0. (Figure taken from Aharo- 
nian etal. [2006a].) 
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with a photon index (F) ranging from 1.7 to 2.9 for sources of different strengths can precisely 
be determined with Monte Carlo simulations with a F of 2.0. Fig. 5.24 shows the ratio of the 
reconstructed to the true A for energy spectra with F ranging from 1.1 to 3.2. Near the safe 
energy threshold at 440 GeV, the differential flux for a source with a F of 2.6 (Crab Nebula) is 
overestimated by 5%, while the differential flux for a F of 1.5 is underestimated by 4%. For 
energies well above the threshold, the energy bias is less than 5% for a wide range of photon 
indices. Thus the A of the reconstructed energy are well suited for spectroscopy of most galactic 
and many extragalactic sources. 

5.4.3 Energy Spectrum of the Crab Nebula 

To verify the precision of the spectral reconstruction, the energy spectrum of the Crab Nebula 
data was analyzed. The results were obtained using the extended cuts, the region-background 
model and the effective area Aext- The differential flux histogram with six bins per decade and 
a fit range from 0.413 GeV to 73.4TeV is shown in Fig 5.25. The fit to a simple power law 
(Eqn. 5.20) is also shown. A good agreement with the results from other experiments is found, 
which are summarized in Tbl. 5.6. Deviations from a power law are apparent at higher energies, 
as expected for an exponential cutoff in the Crab Nebula spectrum, as observed by Aharonian 
et al. [2006a], HEGRA collaboration, R A. Aharonian et al. [2004] and Hillas et al. [1998]. 



Figure 5.25: Energy spectrum of the Crab 
Nebula, fit by a power law. The spectrum 
is determined using the extended cuts, the 
region-background model and the effective 
area Agxt- The exponential cutoff is appar- 
ent at higher energies. 
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Table 5.6: Energy spectra of the Crab Nebula (including statistical errors only), determined from a 
fit to a power law by different experiments: ^ this work, ^ HEGRA collaboration, F. A. Aharonian et 
al. [2004], 3 Hillas et al. [1998] and ^ Masterson and The CAT Collaboration [2001]. 



Experiment 


^ITeV 


F 


4>(£ >lTeV) 




[10""cm-2s-iTeV-i] 




[10-"cm-2s-i] 


H.E.S.S.^ 


2.84±0.05 


2.61±0.02 


1.76±0.04 


HEGRA^ 


2.83±0.04 


2.62±0.02 


1.74±0.04 


Whipple^ 


3.20±0.17 


2.49±0.06 


2.1 ±0.02 


CAT^ 


2.20±0.05 


2.80±0.03 


1.22±0.03 
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Chapter 6 

Imaging with H.E.S.S. 



With the H.E.S.S. data, the first spatially-resolved TeV y-ray sky maps of many y-ray sources 
became available. These maps provided a new picture of many TeV y-ray sources and therefore 
obtained a lot attention (cf. Aharonian et al. [2004b, 2005c]). The H.E.S.S. y-ray maps provide 
images at the highest, previously not attainable energy range, and consequently information 
for the theoretical understanding of y-ray sources. An example is the sky map of RX J1713.7- 
3946 and its relations to the origin of cosmic radiation (Aharonian et al. [2004b]). This chapter 
discusses imaging with H.E.S.S., its limitations and the application of smoothing and image 
deconvolution. 



H.E.S.S. y-ray maps were introduced in the previous chapter. While it is desirable to obtain 
these maps with high resolution, the later is limited by the angular resolution of the shower 
reconstruction and the event statistics. 

6.1.1 Angular Resolution 

H.E.S.S.' angular resolution is limited by the accuracy of shower reconstruction, which deter- 
mines the point spread function (PSF). Mathematically the count map (/) is the result of the 
convolution of the true image (O) with the PSF (PSF), i.e. 



where x and y are the coordinates of the two-dimensional image and * is a short hand notation 
for the convolution operator (Starck et al. [2002]). 

The convolution of a H.E.S.S. sky map can be studied using simulations. Fig. 6.1 shows 

a y-ray map (O) which represents a true y-ray signal and a constant background. The signal 
is modeled with a two-dimensional Gaussian distribution given by Eqn. C.2 with the standard 
deviations = 0.04° and Oi =0.11° in the longitudinal (A) and latitudinal (/3) directions. 
The signal region {0 < 0.3°) contains 4000 y and 14000 background events. The peak intensity 
is 19.3 counts, where 14.3 counts can be attributed to the signal and the rest to the uniform 



6.1 Limitations of H.E.S.S. 7-Ray Maps 




(6.1) 



{P*0){x,y), 



(6.2) 
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background of 5.0 counts/bin. These numbers where chosen according to Tbl. 8.3 to provide a 
realistic simulation. The map has a bin size of 0.01° x 0.01°. 

Fig. 6.2 shows the map / as a result of the convolution of the map O with the PSF. The 
PSF was modeled with the double Gaussian parameterization no. 3 of Tbl. 5.1. The excess 
in / is significantly broader after convolution. This is also reflected by the peak value. After 
background subtraction, the signal peak intensity is only ~6 counts, which is less than half the 
peak intensity of the true excess in O. A profile of the PSF is indicated in the bottom right 
comer. 

6.1.2 Event Statistic 

Another problem of the H.E.S.S. 7-ray maps is their limited event statistic. The total number 
of events is a product of observation time, sensitivity and strength of the source. Sky maps of 
the strongest sources consist of at most a few thousand events. Therefore they are affected by 
Poisson noise in counting statistics. The level of noise in each bin is determined by the number 
of counts. Fig. 6.3 shows a simulation of a count map, which was obtained from Fig. 6.2 by 
replacing each bin content by an integer value given by the Poisson statistic. This count map 
shows the combined effects of convolution and statistical noise. Its difference from the true map 
O is obvious. The intensity profile of the PSF is indicated in the bottom right comer. 

6.2 Image Smoothing 

To produce clearer images Poisson noise is often reduced by smoothing. Smoothing is an effec- 
tive method for reducing the statistical noise through weighted averaging between neighboring 
bins. Here the weights are chosen according to a Gaussian smoothing function, which is a com- 
mon choice. If the smoothed map is 7, the true map O and the smoothing function PSF, then 
smoothing is also expressed by Eqn. 6.2. This equation illustrates that smoothing is a convolu- 
tion. In this work, the integral of PSF is normalized to one so that the total number of counts 
is preserved. A smoothed map appears less noisy and more homogeneous, but also has a lower 
resolution and fewer details because a convolution is implied. 

Fig. 6.4 shows the simulated count map of Fig. 6.3 smoothed with a Gaussian function 
where Gs = 0.03°. The excess appears much more clearly and has a more regular contour. It 
can therefore be identified more easily. A fit of the excess by the Gaussian function G (Eqn. C.6) 
yields the standard deviations of — 0.09° and a; = 0. 14°, which is more than twice the initial 
value of the tme emission in O. 
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Figure 6.1: Simulated map (O) of the true emis- 
sion with a Gaussian signal and constant back- 
ground of 5 counts/bin. 
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Figure 6.2: Map of Fig. 6.1 convolved with the 
PSF. A profile of the PSF is shown at the bottom 
right. 
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Figure 6.3: Simulated count map (/) including 
Poisson noise corresponding to Fig. 6.2. The 
PSF is shown at the bottom right. 
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Figure 6.4: Count map of Fig. 6.3 smoothed 
with a Gaussian function (a, = 0.03°). The PSF 
is shown at the bottom right. 
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6.3 Image Deconvolution 

An alternative method for reducing the noise in sky maps is image deconvolution. It is a tech- 
nique which does not imply a convolution, and is in fact defined as the operation which inverses 
the convolution. Hence, deconvolution appears to be a promising technique for the restoration 
of images which are affected by noise and a limited resolution. Mathematically, image decon- 
volution is defined as solving Eqn. 6.2 for O. The relation between the count map I, the PSF 
(PSF) and the noise term is 

r-OO /"OO 

I{x,y) = / PSF{x-xi,y-yi)0{xi,yi)dxi,dyi+N{x,y) 

Jxi=—°° Jyi=—oo 

= iP*0){x,y)+Nix,y) (6.3) 

(cf. Starck et al. [2002]). Unfortunately it is not possible to determine the exact solution of O in 
the presence of noise and with a limited knowledge of the PSF. However, even the approximate 
solutions can often provide very valuable results. 

6.3.1 The Richardson-Lucy Algorithm 

To find an approximate solution for O in Eqn. 6.3, several methods have been developed. For 
Poisson noise, as in the H.E.S.S. count maps, the Richardson-Lucy (RL) algorithm is used 
and has become a well-established approach. It is a numerical, iterative method for finding an 
approximate solution for O and was first proposed by Richardson [1972] and Lucy [1974]. The 
RL algorithm is defined through the iteration 

0'+\x,y) = 

where / denotes the ith step of iterations and PSF^ the transpose of PSF (Starck et al. [2002]). 
A detailed discussion of the RL algorithm is found in the book by Bertero and Boccacci [1998]. 

The RL algorithm has some important properties which make it very useful. One of these 
properties is the conservation of the total image intensity, which guarantees that the flux at any 
step in the iterations is conserved. Another property guarantees a positive intensity in any region 
of the image and any step of iteration. Therefore unphysical results with negative intensities are 
excluded. As common as it is for iterative deconvolution techniques, the RL algorithm also has 
the property of semi-convergence. Semi-convergence means that there exist an optimal number 
of iterations /opt for which the difference between the restored images 0'°p' and the true image 
O reaches a minimum, /opt depends on individual imaging situations. Therefore, no simple 
stopping rule is known to decide when the optimal number of iterations is reached. However 
/opt can often be determined from simulations. 

6.3.2 Application to H.E.S.S. Data 

To investigate the use of the RL algorithm for the H.E.S.S. data, the RL algorithm was applied 
to H.E.S.S. count maps. The deconvolution of the count maps was done numerically with GDL, 
the "GNU Data Language" (Schellens [2006]), which is a free interpreter of DDL, the "Interac- 
tive Data Language" (Research Systems, Inc. (RSI) [2006]). The RL algorithm was provided 



{PSF*0')(x,y) 



*PSF^{x,y) 



0\x,y), 



(6.4) 
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by the Max_Likelihood procedure from the "IDL Astronomy Library" (Landsman [2004]) 
shown in App. F. This procedure requires convolutions which require Fourier transforms at each 
step of iteration. A Fourier transform is efficiently realized with the "Fast Fourier Transform" 
algorithm. It can be applied if the size of the data set is given by powers of two. Therefore, 128 
bins have been chosen for both dimensions of the count maps. The bin size is 0.01° x 0.01° and 
the scale represents y-ray counts. The necessary PSF parameters were taken from the double 
Gaussian parameterization in Tbl. 5.1. 

Deconvolution of the Sky Map from the Crab Nebula 

As a first example, the RL algorithm was applied to the count map of the Crab Nebula data. The 
analysis of this data was discussed in Chp. 5. The count map is shown in Fig. 6.5 (upper left). 
It consists of ~7000 excess events (Tbl. 5.4) and represents a strong point source. The three 
other count maps show the restorations of the RL algorithm after 10, 40 and 400 iterations. The 
statistical noise of the restored images is clearly reduced. The width of the emission region 
reduces with the number of iterations to a small fraction of its initial width. The count map 
smoothed by convolution with a Gaussian function ((7$ = 0.03°) is shown in the upper right. It 
also shows a clearly reduced noise but with an increased width. The centroid of the excess from 
Tbl. 5.5 is indicated by the black cross in all maps as a reference position. It coincides with 
the centroid of the restored signal. The systematic pointing error of 20" corresponds to half 
the bin width. The H.E.S.S. PSF is indicated in all plots in the bottom right corner. The scale 
represents 7-ray counts. The amplification of noise is negligible and statistical artefacts are not 
observed. The example of the Crab Nebula demonstrates the efficiency of the RL algorithm to 
smooth statistical fluctuations and to mitigate the influence of the PSF in H.E.S.S. count maps 
of strong point sources. 

Deconvolution of the Sky Map from MSH 15-52 

As a second example the RL algorithm was applied to two count maps of data from MSH 
15—52. The analysis of this data is described in Chp. 8. The first count map is obtained with 
the standard image amplitude (M) cut of 80 p.e. and the second with a I A cut of 400 p.e. The 
corresponding energy thresholds are 280 GeV and 900 GeV, respectively. The statistics of the 
two maps are given in Tbl. 8.3. The restoration with the RL algorithm is shown in Fig. 6.6 and 
Fig. 6.7 for 5, 10, 20 and 40 iterations. Again, the statistical noise is significantly reduced. The 
width of the emission region reduces with the number iterations, but the background fluctua- 
tions increase. For comparison, the count map and the smoothed map after convolution with 
a Gaussian function (oi = 0.03°) are also shown in the upper right. The H.E.S.S. PSF is in- 
dicated in all plots in the bottom right comer. The scale represents y-ray counts. The position 
of PSR B 1509— 58 is shown by the black circle. The systematic pointing error of ±0.005° is 
comparable to the bin width. Again the examples of the MSH 15—52 count map demonstrate 
the efficiency of the RL algorithm to smooth statistical fluctuations and to reduce the width 
of an image when applied to H.E.S.S. count maps of strong, extended sources. However, the 
figures with many iterations, especially Fig. 6.7, show a distinctive appearance which could be 
questioned for its physical reality. This question will be addressed in the error analysis of the 
next section. 
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Figure 6.5: Count map of the Crab 
Nebula (point source) restored with the 
Richardson-Lucy algorithm for different 
numbers of iterations (/). With increas- 
ing /, the width of the excess reduces to 
only fraction of its initial value. The orig- 
inal count map and smoothed count map 
are shown for comparison at the upper left 
and upper right. The centroid of the excess 
of the count map is indicated as a reference 
position by the black cross. A profile of the 
PSF is indicated at the bottom right. 
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Figure 6.6: Restored count map of MSH 15—52 (M >80p.e.) with the RL algorithm for different 
numbers of iterations (/). The position of PSR B1509— 58 (black circle) and the PSF are indicated. 
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Figure 6.7: Restored count map of MSH 15-52 (M >400p.e.) with the RL algorithm for different 
numbers of iterations (/). The position of PSR B1509— 58 (black circle) and the PSF are indicated. 
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6.3.3 Error Analysis 

The error of the restoration was determined from forward-folding of simulations. First the true 
map O was simulated to represent emission and background in the source region. This map 
was then convolved with an appropriate model of the PSF. Then noise was added by scattering 
each bin content according to Poisson statistics in order to obtain the simulated count map /. 
These steps are the same as described in Sec. 6.1. The simulated count map was then restored 
using the RL algorithm to generate a restored map. To distinguish between the true signal and 
the statistical artifacts amplified in the restoration, the deconvolution was repeated 20 times for 
each simulated true image O with different random numbers for the noise. 

The true image (O) was then adjusted to minimize the difference between the restored sim- 
ulated count maps and the restored data. The event statistics were provided by the data analysis 
and taken from Tbl. 8.3. In the case of the 80p.e. map, the shape of O was already defined by 
the mathematical model of the fit function. In the case of the 400 p.e. map, O was iteratively 
adjusted until it resembled the restored data sufficiently well. The bin size of the simulated 
maps and the maps of the data was the same. 

The comparison of the 20 different restored maps with each other and with the true map O 
provided an estimate of the error of the restoration. The following quantities are useful for the 
error analysis. 



The Restoration Error 

The restoration error E\x,y) of the pixel {x,y) after i iterations is given by the difference be- 
tween the true and the restored image as E\x,y) — 0{x,y) — j). Since here N — 20 differ- 
ent simulations with different random numbers of noise are compared, 0'{x,y) was replaced by 
the arithmetic mean of the simulations 0\x,y) = 1 /N^^n 0'{x,y) to provide the mean restora- 
tion error E\x,y). E\x,y) is then given as 

E\x,y) = 0{x,y)-0\x,y). (6.5) 

E\x,y) expresses the mean difference between the true and the restored map. 



The Standard Deviation of the Restoration 

Similar to the restoration error, the standard deviation of the restoration for each pixel S'{x,y) 
can be defined. If N is the number of simulations with different random numbers, then 

S-(,,).(fI^fc5W)l!)\ (6.6) 

S\x,y) represents the noise of each pixel of the restored map. 



Distribution of Noise of the Restoration 

The noise in the restoration is determined by the deviation from the mean restoration value. The 
deviation of each pixel D\x,y) from the mean restoration value is given as 

D\x,y) = 0\x,y)-0\x,y). (6.7) 
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The distribution of Z)' for all pixels {x,y) is a measure of the overall noise distribution of the 
restoration at step In contrast to S\x,y), the deviation does not refer to a particular pixel or 
region but to the deviations within the whole map. The distribution is a measure of the accuracy 
of the restoration. 



The Relative Norm of the Restoration Error 

According to [Bertero and Boccacci, 1998, pg. 110], the relative norm of the restoration error 
(g/) is defined as 

\\0(x,y)\\ 

where | |0(jc, 1 1 is the Euclidean norm of 0{x,y) defined as 



\\0{x,y)\\=mO{xk,yi)n . (6.9) 
\k,i J 

Hence e, is a positive number which expresses the difference between the true (O) and the 
restored image {O') after i iterations. In the case of a perfect restoration, e, = 0. Due to the 
semi-convergence of the RL algorithm, the sequence of £, has a minimum. It defines the opti- 
mal number of iterations (iopt) which provides the smallest relative restoration error (Copt). In 
principle, is an ideal measure for the agreement between the true image and the restoration. 
However, for images with compact signal regions and non-negligible background, £, varies with 
the size of the image. The reason for this is that a variation of the image size will only change 
the size of the background area, and the deconvolution of background will only amplify noise 
and hence increase the restoration error. Therefore, other measures such as those discussed 
above are also necessary to assess the restoration. 

Since N = 20 simulations with different random numbers of noise are considered here, is 
used as the arithmetic mean of the individual 



Errors of the 80p.e. Model 

The true map O of the 80p.e. map is modeled with a Gaussian distribution (Gw = 0.04°, O/ = 
0. 1 1°) according to the fit of the corresponding excess map (Tbl. 8.4). The statistics in the signal 
region with a radius 6 of 0.3° is given by 4000 excess and 14000 background events (Tbl. 8.3), 
yielding a peak excess of 14.3. The simulation is the same as the one discussed in reference to 
Fig. 6. 1 and was repeated 20 times with different random numbers of statistical noise. The PSF 
was chosen according to parameterization no. 3 of Tbl. 5.1. Since the performance of the RL 
algorithm is invariant under rotations of an image, the orientation of the excess distribution was 
not adjusted in the simulations. The width and length of the excess are oriented along the j8 — 
and A —axis, respectively. 

The restoration of one of the 20 simulated count maps is shown in Fig. 6.8 for 5, 10, 20 and 
40 iterations. It is similar to the restored count map of the data in Fig. 6.6. The width of the 
emission region reduces similarly. For comparison, the last plot shows the count map that was 
smoothed by convolution with a Gaussian function (Os = 0.03°). It has the largest width. 

Fig. 6.9 shows the mean error E\x,y) (Eqn. 6.5) of the 20 different restorations. The scale is 
chosen in percentages of the peak intensity of the simulated excess of 14.3 counts. The largest 
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error {E'^^^) is found at the center of the excess, where the restored intensity is systematically 
too small. On both sides of the center region, at a distance of about 0.1° to 0.2°, the excess is 
systematically too high. A very similar error distribution is seen in the smoothed map. E^x^y) 
decreases with an increasing number of iterations and the restoration approaches the true image 
O. For example, after 40 iterations E^^^-^ is smaller than ~15%. In comparison, the E',^ of the 
smoothed count map shows the largest error of ~65%. 

Fig. 6.10 shows the standard deviation S\x,y) (Eqn. 6.6) as obtained from the 20 different 
simulations. Again, the scale is chosen in percentages of the peak amplitude of the simulated 
true emission. An increase of the standard deviation with the number of iterations can be seen. 
The largest deviations are found at the center region. After 40 iterations they reach a level 
of ~10%. For comparison, the standard deviation in the smoothed count map is also shown. It 
has the smallest deviations — only ~3%. 

Fig. 6.11 shows the distributions of the pixel's noise D' {x, y) (Eqn. 6.7) for different numbers 
of iterations /. The distributions of the 20 simulations are added to one plot. The noise is 
measured in percentages of the peak intensity of the excess of the true image O. The noise of 
the restoration is approximately Gaussian distributed with a zero-mean as shown by a fit to a 
Gaussian function (red line). Only for higher numbers of iterations are small deviations from a 
Gaussian distribution found. The level of noise can be characterized by the standard deviations 
of 1.3, 2.2, 3.4 and 5.3% for 5, 10, 20 and 40 iterations, respectively. The noise of the map after 
5 iterations is comparable to the noise of the smoothed map. 

The last plot of Fig. 6.11 shows the mean relative norm of the restoration error (£, ) (Eqn. 6.8). 
Ei is calculated for the full region of 128 x 128 bins and as the mean of the 20 different sim- 
ulations of the 80p.e. map (red marker; the black marker refers to the 400 p.e. map which 
is discussed later). The minimum is obtained after about 10 iterations. The statistical error is 
~2%. Within the range of about ±5 iterations from the minimum, Copt changes by less than 
5%, showing that the restorations change only marginally in the interval 5 < / < 15. 

Fig. 6.12 shows the corresponding excess profiles of the restored maps of Fig. 6.8 after 
background subtraction. The distributions show the width (black) and length (blue) along the 
A- and j8-axis of the excess, respectively. The solid line represents the Gaussian fit function 
(Eqn. C.2). The fit value of the standard deviations (a^, a/) of the Gaussian function is shown 
in each plot. For example, for 5, 10, 20 and 40 iterations, the width represents 180, 155, 130 and 
120% of the simulated excess distribution in O, while the smoothed map (plot 5) has a width of 
220%. Small waves in the profiles show the amplification of noise increasing with the number 
of iterations. The smoothed map shows only little noise. The last plot shows the successive 
decrease of the width and length versus the number of iterations. 

The individual errors discussed here are summarized in Tbl. 6.1. 
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Figure 6.8: Restored simulated count 
map of Fig. 6.1 (M >80p.e.) with the 
Richardson-Lucy algorithm for different 
numbers of iterations (/). A profile of the 
PSF is indicated at the bottom right. With 
increasing /, the width of the emission re- 
gion reduces and a slight asymmetry of the 
excess appears. The smoothed count map 
is shown for comparison. 
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Figure 6.9: The mean restoration error 
{E'{x,y)) between the restored and the 
simulated image (O). The mean is ob- 
tained from the restoration of 20 different 
count maps. The scale is in percentages 
of the peak intensity of the simulated ex- 
cess. E'{x,y) decreases with increasing it- 
erations, while background fluctuations in- 
crease. The largest deviations are found in 
the center region. E'{x,y) is highest for the 
smoothed count map, which is shown for 
comparison. 
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Figure 6.10: Standard deviation of the 
restoration {S'{x,y)) of 20 different re- 
stored count maps. The scale is given in 
percentages of the peak intensity of the 
simulated excess. S'{x,y) increases with 
the number of iterations. The highest val- 
ues are reached in the center region. Af- 
ter 40 iterations, S'{x,y) reaches a level of 
~10%. S'{x,y) of the smoothed count map 
is comparatively small. 
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Figure 6.11: Distribution showing each pixel's deviation from its expectation value. The deviation 
is measured in percentages of the peak intensity of the simulated excess and represents the noise. 
A fit to a Gaussian distribution is shown (red). The last plot shows the mean relative norm of the 
restoration error (£,) versus the number of iterations /. The minimum is found for /opt ~ 10 iterations. 
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Figure 6.12: Profiles of the restored count maps of Fig. 6.8 showing the width and length of the 
excess. With increasing number of iterations (/), the width and length decreases while the amplifica- 
tion of noise increases. The width and length of the smoothed count map are shown for comparison. 
The bottom right plot shows the width and length measured in standard deviation vs. the number of 
iterations. 
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Restoration of the 400 p.e. Map 

The curved excess of the 400 p.e. map was modeled by an arc consisting of all points on a 
ring (0.14° < r < 0.16°) around the origin which lie in the fourth quadrant of the coordinate 
system. In addition, the arc is convolved with a Gaussian function with a standard deviation of 
0.02° and shifted to the middle of the map. This model was found through a comparison of the 
restored count maps of the simulation and those of the data. The statistics in the signal region 
with a radius 6 of 0.3° are given by 1000 excess and 400 background events (Tbl. 8.3), yielding 
a peak excess of 7 counts in O and a background of 0.14 counts/bin. The simulated maps are 
shown in Fig. 6.13 to 6.16. The count map is only one of 20 with different random numbers of 
noise which have been restored. The PSF was modeled according to parameterization no. 5 in 
Tbl. 5.1. 

A typical restoration of one of the 20 simulated count maps with the RL algorithm is shown 
by the images in Fig. 6.17 for 5, 10, 20 and 40 iterations. The restorations are similar to 
the restoration of the count map of the data in Fig. 6.7. The width of the emission region 
reduces accordingly. For comparison, the last plot shows the count map that was smoothed by 
convolution with a Gaussian function (oi = 0.03°). It shows the widest extension. 

Fig. 6.18 shows the mean error E\x,y) (Eqn. 6.5) of the 20 different restorations. The scale 
is chosen in percentages of the peak intensity of the simulated excess (7 counts). The largest 
error (E^^^) is found along the arc, where the restored intensity is systematically too small. 
On both sides close to the arc, the restored excess is systematically too high. A very similar 
error distribution is seen for the smoothed map. E^ (jc, y) decreases with an increasing number of 
iterations as the restoration approaches the true image O. For example, after 40 iterations, E'^^ 
is smaller than ~15%. In comparison, the E]^^ of the smoothed count map shows an error of 
~65%. 

Fig. 6.19 shows the standard deviation S\x,y) (Eqn. 6.6) as obtained from the 20 different 
simulations. Again, the scale is chosen in percentages of the peak amplitude of the simulated 
true emission. An increase of the standard deviation with the number of iterations can be seen. 
The largest deviations are found at the arc. After 40 iterations, they reach a level of ~30%. 
For comparison, the standard deviation in the smoothed count map is also shown. It has the 
smallest deviations — only ~3%. 

Fig. 6.20 shows the distributions of the pixel noise D\x,y) (Eqn. 6.7) for different numbers 
of iterations /. The distributions of the 20 simulations are added to one plot. The noise is 
measured in percentages of the peak intensity of the excess of the true image O. The distribution 
is fitted with by a Gaussian function with a zero-mean (red line). The noise is clearly not 
Gaussian distributed. The standard deviations of the distributions are 1.1, 1.5, 2.2 and 3.1% for 
5, 10, 20 and 40 iterations, respectively. The smoothed map has the smallest standard deviation 
of 0.7%. The deviations from a Gaussian distribution of noise is an indication of statistical 
artefacts which are created in the restoration. 

The last plot of Fig. 6.11 shows the mean relative norm of the restoration error (e,) (Eqn. 6.8). 
Si is calculated for the full region of 128 x 128 bins and as the mean of the 20 different simu- 
lations of the 80 p.e. map (black marker). The minimum is obtained after about 10 iterations. 
The statistical error is ~2%. Within the range of about ±5 iterations from the minimum, Copt 
changes by less than 5%, showing that the restorations for 5 < / < 15 will provide restoration 
errors similar to £opt. 

The individual errors discussed here are summarized in Tbl. 6.1. 
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Figure 6.13: Simulated map (O) of the true 
emission with arc-like signal and constant back- 
ground of 0.14 counts/bin (M >400p.e.). 
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Figure 6.14: Map of Fig. 6.13 convolved with 
the PSF. A profile of the PSF is shown at the 
bottom right. 
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Figure 6.15: Simulated count map (/) including 
Poisson noise corresponding to Fig. 6.14. The 
PSF is shown at the bottom right. 
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Figure 6.16: Count map of Fig. 6.15 smoothed 
with a Gaussian function (a, = 0.03°). The PSF 
is shown at the bottom right. 
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Figure 6.17: Restored simulated count 
map of Fig. 6.13 (M >400p.e.) using the 
Richardson-Lucy algorithm for different 
numbers of iterations (/). A profile of the 
PSF is indicated at the bottom right. With 
increasing /, the width of the emission re- 
gion reduces and the noise increases. The 
smoothed count map is shown for compar- 
ison. 
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Figure 6.18: Mean restoration error 
{E'{x,y)) between the restored and the 
simulated image (O). The mean is ob- 
tained from the restoration of 20 different 
count maps. The scale is in percentages 
of the peak intensity of the simulated ex- 
cess. E'{x,y) decreases with increasing it- 
erations; background fluctuations increase. 
E'{x,y) is greatest for the smoothed count 
map, which is shown for comparison. 
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Figure 6.19: Standard deviation {S'{x,y)) 
of 20 different restored count maps. The 
scale is given in percentages of the peak 
intensity of the simulated excess. With an 
increasing number of iterations, S'{x,y) in- 
creases. The highest values are reached 
in the center region. After 40 iterations, 
S'{x,y) reaches a level of ~30%. S'{x,y) of 
the smoothed count map is comparatively 
small. 
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Figure 6.20: Distribution showing each pixel's deviation from its expectation value after certain 
numbers of iterations /. The deviation is measured in percentages of the peak intensity of the simu- 
lated excess and represents the noise. A fit to a Gaussian distribution is shown (red). The last plot 
shows the mean relative norm of the restoration error (e,) versus the number of iterations /. The 
minimum is found at j'opt ~ 10 iterations. 
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6.4 Conclusion 



Table 6.1: Errors of the restoration of the count map of MSH 15—52 (M >80p.e.). The values E'^^, 
E'^ and 00,1 are measured in percentages of the peak excess counts. refers to the width of the 
excess of the true map O. 
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Table 6.2: Errors of the restoration of the count map of MSH 15—52 (lA >400p.e.). The values 
^max' "^max ^D,i are measured in percentages of the peak excess counts. 
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Systematic Errors of the Simulations 

While previous discussion refers to the statistical errors produced by Poisson noise, it is also 
possible to estimate the systematic errors of the restored sky maps. The major source for sys- 
tematic errors has been found in deviations of the morphology. The influence of the morphology 
on the restorations has been investigated with simulations. Before, unphysical images were ex- 
cluded and the set of possible maps of O were reduced to a realistic subset of smooth maps 
with simple geometric structures. With these restrictions, only small variations were possible in 
order to reproduce count maps of similar appearance. The parameters which were investigated 
included the size as well as the aspect ratio of the excess distributions. The results suggest that 
an estimate for the systematic error would be less than ~10% and ~20% of the peak intensity 
of the excess for the 80 and 400 p.e. maps, respectively. 

Other potential sources of systematic errors have also been investigated, including the num- 
ber of excess and background events and the size of the PSF. No systematic errors related to 
these parameters have been found. 



6.4 Conclusion 

Two methods of image restoration for H.E.S.S. count maps have been discussed here: smooth- 
ing by convolution with a Gaussian function and the Richardson-Lucy algorithm for image 
deconvolution. Both methods can significantly reduce the statistical noise in count maps and 
reveal morphological details which are hidden by statistical noise. 
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Convolution is a straightforward approach which provides stable results which are only 
little affected by statistical fluctuations. Simulations show that the smoothed map for the 80 and 
400p.e. maps of MSH 15—52 have a restoration error of about (65±3)%. 

On the other hand, the Richardson-Lucy algorithm is more sensitive to noise in count maps 
and therefore requires simulations for error analysis. With simulations, it is possible to estimate 
the quality and error of the restorations and to limit the statistical fluctuation to an acceptable 
level. At the expense of increasing noise, the RL algorithm can provide a high restoration of 
morphological details. Depending on the objective of the analysis, one can choose between 
a restoration with less details and small errors or with higher details but also more statistical 
artefacts. For the analysis of the morphology of MSH 15—52 in this work, preference was given 
to smaller errors. 10 iterations with the RL algorithm were considered a good compromise. In 
the case of the 80p.e., this provides a stable restoration with a mean error of at most (30±5)% 
from the true value. The noise in the restored map is Gaussian distributed and has a standard 
deviation of 2.2%. The 400 p.e. map after 10 iterations with the RL algorithm provides a 
relative error of the excess of at most (30±13)%. These maps are in agreement with the maps 
obtained by smoothing. The simulations of the y-ray maps of MSH 15—52 have shown that the 
application of the RL can be useful and provide y-ray maps of high resolution. 
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Chapter 7 

Search for Pulsed Emission from Pulsar 



When 7 radiation from pulsars is detected, the question arises whether a part of this radiation 
is pulsed. The detection of pulsed 7 radiation from pulsars is of special interest for the de- 
velopment of pulsar models, since it can provide information about the emission process and 
constrain parameters. Unfortunately the detection of pulsed TeV 7 radiation from a pulsar wind 
nebula is more difficult since the nebula often emits a significant amount of constant 7 radiation. 
If at all, the pulsed flux will only be observable as one component of the 7-ray flux from the 
source region. This chapter will introduce methods for the analysis of H.E.S.S. data for pulsed 
emission, statistical tests for evaluation and the calculation of flux upper limits. 

7.1 Pulsar Light Curves and Phasograms 

Methods for finding pulsed emission from pulsars rely on an analysis of the pulsar light curve. 
The pulsar light curve can be represented in a phasogram. A phasogram shows the light curve 
during one period, usually obtained as an average over many periods. To produce a phasogram 
from H.E.S.S. data, the standard analysis is first carried out to reconstruct the 7-ray showers as 
described in Sec. 5. Then the events contained in a small signal region which encompasses the 
pulsar are filled in an phasogram. Since the event statistic is limited, between 16 to 25 bins are 
chosen. The phase of each event is calculated from its time stamp and the pulsar ephemeris. A 
H.E.S.S. phasogram for the Crab Pulsar is shown in Fig. 7.2. 

7.2 Ephemerides 

Pulsar ephemeris contain a set of parameters which can describe the pulsar phase as a function 
of time Since pulsars behave like a slowly decelerating sphere, the first terms of the 

Taylor series 

(0 = 00 + / • - To) + ^f-{t-Tof (7.1) 

already provide a description of the pulsar phase with high precision. The parameters in the 
Taylor series already constitutes an ephemeris. In this equation ^0 is the phase at the reference 
time Tq and / is the frequency of rotation of the pulsar. Since the pulsar phase is often very 
precisely determined by radio observation, ephemeris are usually derived from radio measure- 
ments. Given the ephemeris of a pulsar, one can calculate its phase at any time (t). However, 
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pulsars do not exactly obey Eqn. 7.1 — they are also subject to other processes such as glitches. 
Therefore, an ephemeris is only valid within a certain validity period which can range from a 
few days, as in the case of the Vela Pulsar, to many years, as in the case of PSR B 1509— 58. 
The validity period is therefore also part of the ephemeris. There are some other parameters 
necessary for a precise definition of an ephemeris. An example is shown in Tbl. 7.1. 

Several formats for ephemeris exist. A common format for pulsar ephemeris is the GRO 
format. It was first used by the Compton Gamma-Ray Observatory community and is now the 
standard format of the Australian Pulsar Timing Data Archive (Manchester et al. [2006]) from 
the Australia Telescope National Facility (ATNF). Since several ephemerides have been taken 
from this archive for the analysis of H.E.S.S. data, this format was implemented to the H.E.S.S. 
standard analysis. The GRO format consists of one line containing 13 parameters in the case of 
a single pulsar and of two lines with 10 additional parameters in the case of a binary system. 
An illustration of the ephemeris format and parameters is given in Tbl. 7.1. 



Table 7.1: Structure and contents of the GRO ephemeris format as documented at Manchester et al. 
[2006]. The basic parameters are required for any pulsar. The binary parameters are additional 

parameters for pulsars in binary systems. 



Character 


Symbol 


Meaning 


Basic Parameters 


1-8 




Pulsar name (truncated if a J2000 name) 


10-21 


a 


J2000 right ascension [hh mm ss.sss] 


23-34 


5 


J2000 declination [-dd mm ss.ss] 


36-40 




Start of vaUdity range [MJD] 


42-46 


Tmax 


End of validity range [MJD] 


48-62 


tGEO 


TDB epoch of pulse frequencies and infinite frequency UTC 






pulse TOA at geocenter [MJD] 


64-80 


f 


Pulse frequency at the Solar system barycenter [Hz] 


82-93 


f 


1st time derivative of bary centric pulse frequency [s~^] 


96-104 


f 


2nd time derivative of barycentric pulse frequency [s~^] 


106-109 




RMS residual of fit in milliperiods 


111 




Letter code indicating origin of data (A = Australia) 


115-119 




Planetary system ephemeris used for barycenter correction 


121-130 




Full J2000 pulsar name 


Binary Parameters 


1-8 




Pulsar name (truncated if a J2000 name) 


10-25 


Pb 


Orbital period (at the Solar system barycenter) [s] 


26-37 


X 


Semi-major axis of pulsar orbit [s] 


39-48 


e 


Orbital eccentricity 


50-63 


To 


TDB epoch of periastron passage [MJD] 


65-74 


CD 


Longitude of periastron [deg] 


76-82 


6) 


Rate of periastron advance [deg/yr] 


84-91 


r 


Time dilation and gravitational redshift term [s] 


93-102 


Pb 


First time derivative of orbital period 


104 




Letter code indicating origin of data (A = Australia) 
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7.3 Time of Flight Corrections 

Although the pulsar phase is described by the ephemeris, one cannot directly apply Eqn. 7.1 to 
the event time provided by the time stamp in order to calculate the correct pulsar phase {^{t)). 
The reasons for this are mainly the changing distance between the pulsar and the observatory, 
but also e.g. relativistic effects. Therefore it is necessary to correct the event time before 
Eqn. 7.1 is applied. These corrections are described below. 



7.3.1 Solar System Barycenter Correction 

The solar system barycenter (SSB) correction compensates for the changing distance between 
the pulsar and the observatory caused by the motion of the earth. The situations is illustrated in 
Fig. 7.1. Eqn. 7.1 is correct, if applied to the arrival times at the SSB, since the SSB provides 

a reference frame of constant velocity in good approximation. Arrival times of the same pulse 
recorded at the observatory have a time difference (A?ssb) with respect to the SSB depending on 
the earth's position. The relation between the arrival time at the SSB (tt,) and at the observatory 
(t) is given as 

„ = ,-A,ssB=»-^55^^^, (7.2) 

C 

where is a vector to the SSB from the phase center of the observer, fipsR is a unit vector 
pointing from the SSB to the pulsar and c is the speed of light, is calculated as 

^SSB(0=4(0+?r(0, (7.3) 

where (t) is the vector from the geocenter to the SSB and er{t) is the vector from the geocenter 
to the observer. er{t) is determined by 

(cosA(?)cos<^A 
sinA(?)cos0^ , (7.4) 
sin<^^ / 

where is the earth's radius and (A, 6^) are the geographic coordinates of the observer. The 
coordinates of the SSB and the geocenter are obtained from the DE200 solar system ephemeris 
(Standish [1982]). The DE200 ephemerides provide the earth's position for times between 
the years 1980 and 2020 with an accuracy of within a few meters corresponding to a time 
resolution of nanoseconds, which is more than sufficient for time of flight corrections. Further 
details regarding the coordinate and time transformations for H.E.S.S. can be found in Gillessen 
[2004]. 

The magnitude of (A?ssb) can be estimated with a simple example as follows. A geostation- 
ary observatory is orbiting the sun with an approximate velocity Vearth ~ 30km/s. If a pulsar 
with a typical frequency of about /p = 30 Hz emits pulses, then the distance between two pulses 
is Ap = c//p ~ 10000km. The earth travels this distance in Tp = Ap/vearth = 333 s. Conse- 
quently, if the earth is moving directly towards or away from the pulsar, the pulsar phase can 
shift by one period in about 5 minutes. 



7.3.2 Binary Correction 

If the pulsar is part of a binary system, then phase shifts also arise due to the orbital movement 
of the pulsar around its companion. The phase shift can be measured and the orbit can be param- 
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Source 




Figure 7.1: Illustration of the solar sys- 
tem barycenter correction. Planes of equal 
phases are perpendicular to the line of sight 
of the pulsar. Thus the observed pulse 
phase changes with the position of the 
earth. A similar phase shift is introduced 
by the orbit of a pulsar in a binary sys- 
tem. (Figure taken from [Schmidt, 2005, 
pg. 44].) 



eterized and summarized in the pulsar ephemeris. A mathematical model (BT), which describes 
the time shift as observed at the SSB, was derived by Blandford and Teukolsky (Blandford and 
Teukolsky [1976]). The BT model defines the transformation from the time of arrival at the 
SSB (th) to the pulsar proper time (T). The BT formula contains Keplerian elements common 
for the description of binary systems and relativistic terms as follows: 

tb-to = r + {jcsin(D(cos£' -e) + [xcos(d(1 -e^)^/^ + 7] sinE} X 
{1 [.)ccos(d(1 — e^y^^ cosE — xs'mcosinE] x 

Pb 

[l-ecosE)-^}. (7.5) 

Here Pt,, e and CO are the binary orbital period, orbital eccentricity and longitude of the periastron 
respectively. The longitude of the periastron is defined as the angle between the periastron and 
the ascending node^\ x is the projected semi-major axis^^ of the pulsar orbit in time units. 
7 measures the combined effect of gravitational redshift and time dilation, to is an arbitrary 
reference time. The eccentric anomaly E is defined by Kepler's equation, 

E-esinE = —{th-To), (7.6) 

in which Tq is a reference time of periastron passage, measured in the TDB system. 



7.3.3 TEMPO and CRASH 

The timing corrections for the H.E.S.S. data are done with the Coordinate Transformation Soft- 
ware for H.E.S.S. (CRASH) (H.E.S.S. collaboration [2001]), which is part of the H.E.S.S. anal- 
ysis software. CRASH implements the DE200 solar system ephemeris (Standish [1982]) and 

'^The ascending node is the point in the orbit of an object when it crosses the ecliptic (i.e. celestial equator) 
while moving from south to north. 

^'The projected semi-major axis is the semi-major axis of the apparent ellipse, i.e. of the projection of the actual 
elliptical orbit onto a plane perpendicular to the line of sight of the observer. 
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provides the CRASH pulsar class, which can calculate the pulsar phase at the SSB based on 
pulsar ephemerides. CRASH can read pulsar ephemerides of different formats. The implemen- 
tations were verified with TEMPO (Taylor et al. [2004]), the standard program for pulsar timing 
in radio astronomy. TEMPO can deduce pulsar rotation as well as astrometric and binary pa- 
rameters by fitting models to pulsar data. It also contains an implementation of the BT model. 
Information about the implementation of the GRO format in CRASH and its verification can 
be found in Breitling et al. [2004]. The calculation of the pulsar phase including the SSB cor- 
rection is done in the H.E.S.S. standard analysis chain, which can also produce phasograms of 
arbitrary energy ranges. 

The implemented timing capabilities of the H.E.S.S. standard analysis were tested with 
optical H.E.S.S. data from the Crab Pulsar (PSRB0531+21). In this test, one telescope was 
equipped with a special PMT camera which was developed for the detection of optical light 
curves from the Crab Pulsar. Further details about this experiment can be found in Hinton et al. 
[2006] and Masterson and H. E. S. S. Collaboration [2003]. The ephemeris for this analysis was 
taken from the Australian Pulsar Timing Data Archive and is shown in Tbl. 7.2. Fig. 7.2 shows 
the phasogram with the mean ADC counts of 100 s (2.9 x 10^ events) of optical data taken on 
November 23, 2003 (21:37 UTC). The histogram of 25 bins (grey) is overlaid with a histogram 
with 100 bins which can resolve the pulse shape more clearly. For comparison, the optical light 
curve as obtained by OPTIMA in January 2002 is shown in Fig. 7.3. OPTIMA is a high-speed 
photo-polarimeter for optical pulsar measurements of high precision (Kanbach et al. [2005]). A 
good agreement between the two light curves is found. 

Table 7.2: GRO ephemeris of the Crab pulsar taken from the ATNF archive (Manchester et al. 
[2006]). The parameters are documented in Tbl. 7.1. The validity range is given in Modified Julian 
Dates and corresponds to a range from January 31, 2004 to July 3, 2005. 



Character 


Value 


a 


05 34 31.972 


5 


22 00 52.07 


Tram 


52944 MJD 


T 

-'max 


52975 MJD 


tCEO 


52960.000000296 MJD 


f 


29.8003951530036 s-i 


f 


-3.73414D-10S-2 


f 


1.18D-20S-3 
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Crab Pulsar, Optical Observation by H.E.S.S. 



X^ = 2.55e+06, Prob. <10 
Z§ = 9.44e+06, Prob. <10 " 
H(2) =7.916-29, Prob. <10''' 
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Figure 7.2: Phasogram for 100 s of the opti- 
cal H.E.S.S. data from the Crab Pulsar obtained 
with the standard analysis. The same data is rep- 
resented with a binning of 25 and 100 bins. 




0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.: 
Phase 

Figure 7.3: Optical light curve of the Crab 
Pulsar as measured with high precision by 
OPTIMA. (Figure taken from Kanbach et al. 
[2005].) 



7.4 Statistical Test 

For many pulsar observations with H.E.S.S., in particular for observation of PWNs, the detec- 
tion of pulsed emission requires the identification of a small pulsed component from a strong 
constant background signal. The situation is complicated by the fact that each pulsar has its 
own characteristic light curve which can differ with wavelength (cf. Fig. 2.8). Thus the position 
and shape of the pulse are not known and its identification is more difficult. For detecting an 
unknown light curve above a strong background, statistical tests can be applied. Common tests 
are the x^, the and the H test. They are tests of uniformity of a phasogram. For weak pulse 
shapes the tests differ in their power. The power of a test is its efficiency to reject the hypothesis 
of a uniform phase distribution, if this is not true. A detailed discussion of these tests and their 
application can be found in de Jager et al. [1989] and de Jager [1994]. Here they are briefly 
introduced. 



7.4.1 The Test 

Pearson's test is applied by fitting a constant to the phasogram. The number of degrees 
of freedom is « — 1, where n is the number of bins in the phasogram. The typical choice for 
H.E.S.S. data is 16 < n < 25. The x^ value or rather the corresponding x^ probability then 
determine whether the hypothesis of a uniform distribution can be rejected or not. The X^ test 
is most efficient for single and narrow pulse shapes but less for other, especially sinusoidal, 
shapes. Another disadvantage is the dependence on the binning. 



7.4.2 The Zf^ Test 



The Z^, test can be considered as a complementary test to the x^ test, in the sense that it is 
more sensitive to sinusoidal signals with the periodicity m. Z^ is defined through the sum of 
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trigonometric moments aj and Pj as 



zi = 2NY^{aj + p}), (7.7) 



with 



i£cosa(/),) (7.8) 

!=1 



1 ^ 



and 

Here A'' is the number of events in the phasogram and is the phase of event /. is 
distributed for 2m degrees of freedom. The probability is calculated accordingly, m = 2 is a 
good choice for the detection of wider pulse shapes, which is often used for the analysis of 
H.E.S.S. data. Another advantage of the test, besides sensitivity to sinusoidal signals, is its 
independence from the binning of the phasogram. 



7.4.3 The H test 

The H test is an improved version of the Z^ test which has an increased sensitivity for arbitrary 
pulse shapes. It is very useful when no a priori information about the pulse shape is available. 
The H test is as powerful as the test and more powerful than the Z^ test in the case of more 
than two pulse peaks (cf. de Jager et al. [1989] and de Jager [1994]). H is defined through the 
Z^ test as follows: 

H= max (Zl-4m + 4) (7.10) 

1<OT<20 

The selection of the maximum value of H accounts for its increased power in comparison to the 
Z^ test. In the case of the absence of a signal, H is distributed by an exponential function and 
the probability P for obtaining a lager value for H is determined as 

P{>H) = exp{-0.4H). (7.11) 

This test is a favorable choice for the analysis of most H.E.S.S. pulsar data, since often no priori 
information about the pulse shape is available. 



7.4.4 Application of Tests to the Optical Crab Pulsar Data 

These three tests have been applied to the optical H.E.S.S. data from the Crab Pulsar shown in 
Fig. 7.2. The results are listed in the panel on the figure. For all three tests, the probabilities for 
a uniform light curve are close to zero as is expected for data with such a clear pulse shape as in 
the Crab Pulsar. The numerical values of the probabilities are less than the numerical accuracy 
of 10-19. 
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7.5 Calculation of Upper Limits 

If pulsed emission cannot be detected, one can still determine an upper limit of the pulsed flux 
in a phase region. The phase region is referred to as an On-region and the remainder as an 
Off-region. The On-region is usually chosen according to the light curve at other wavelengths 
where a pulsed signal has been observed. 

Eqn. 5.24 and Eqn. 5.25 can provide the flux 4> and the corresponding error o$ within a 
given energy range. A common choice for H.E.S.S. data is the energy range above the thresh- 
old energy or above 1 TeV. 4> and <7<i> are sufficient to calculate the upper limit of the flux 
(UL^) for a confidence level (CL). There have been different methods proposed for the calcu- 
lation of upper limits. Here the unified approach by Feldman and Cousins [1998] was chosen, 
which is described in more detail in App. D. UL^ is determined with an upper limit function 
Ful{^, <7$,CL) as described in App. D. Examples for the application of the H.E.S.S. standard 
analysis and the calculation of upper limits according to Feldman and Cousins are found in 
Schmidt et al. [2005], where H.E.S.S. data of young pulsars are analyzed. 
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After the introduction and verification of the H.E.S.S. data analysis techniques, they are now 
applied to the H.E.S.S. data from MSH 15—52. The observation, the data, the analysis and the 
results are subject of this chapter. The detection, position, energy spectrum and morphology of 
MSH 15—52 as well as the analysis for pulsed emission from PSR B 1509— 58 are discussed. 



8.1 Observation 

MSH 15—52 was repeatedly observed in 2004 from March 26 to July 20. A total of 78 stereo- 
scopic observation runs with a duration of mostly 28 minutes each were taken during this period. 
The data was taken in wobble mode with all four telescopes fully operational. The system trig- 
ger required a telescope multiplicity of at least two. The observation period was favored by good 
weather providing data of high quality. Fig. 8.1 shows the corresponding count map of y-ray 
candidates after run selection and cuts. The galactic plane is indicated by the yellow line in the 
upper right region. It passes the field of view at the approximate distance of 1° from source 
position. The position of PSR B1509— 58 as determined from radio observations is represented 
by the black circle near the center of the map. 



a 
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a 
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Figure 8.1: Count map of the H.E.S.S. 

data from MSH 15—52 passing the quality 
criteria. White crosses represent the target 
position (center) and the four observation 
positions for the different wobble offsets. 
The position of PSR B1509-58 is marked 
(black circle). The galactic plane is indi- 
cated by the yellow line in the upper right. 
The PSF is shown in the lower right comer. 
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8.1.1 Observation Position 

The data was taken at the target position of (J2000) (15*^14^27% -59°16'18'0 and wobble off- 
sets of ±0.5° in RA and ±0.5° in Dec, resulting in four different observation positions. The 
wobble positions and the target position are indicated by white crosses in Fig. 8.1. Tbl. 8.1 
shows the number of runs after run selection for each wobble position. The number of pos- 
itive and negative wobble offsets is approximately equal for each wobble direction reducing 
systematic errors. The full run list is found in App. E. 

Table 8.1: Summary of observation runs taken at the four different wobble offsets shown in Fig. 8.1. 



Obs. Pos. 


No. of runs 


<Zenith angle> [°] 


live-time [h] 


Dec +0.500° 


23 


36.7 


9.93 


Dec -0.500° 


24 


37.7 


9.83 


RA +0.978° 


7 


36.7 


2.98 


RA -0.978° 


8 


36.7 


3.40 




62 


37.1 


26.14 



8.1.2 Run Selection 

The runs were selected according to the quality criteria described in Sec. 5.3.1. The result of 
the run selection is summarized in Tbl. 8.2. Run by run statistics are given in App. E. The 
low rejection ratio of 16 out of 78 runs (~20%) is owed to good weather conditions during the 
observation period. The total observation time was 35 hours of which 29 hours had sufficient 
data quality. After dead-time correction, data with a live-time of 26. 14 hours remained. 

Table 8.2: Results of the run selection according to the quality criteria of Sec. 5.3.1. 





No. of runs 


Observation time [h] 


selected 


62 


29.1 


rejected 


16 


6.2 


total 


78 


35.3 



8.1.3 Zenith Angle Distribution 

The mean zenith angle of the selected data was 37.1°. Fig. 8.2 shows the zenith angle dis- 
tribution of the 7-ray candidate events. The three peaks result from different zenith angles of 
culmination for each month of the observation period. 

8.2 Detection of the 7-Ray Signal 

The selected data has been processed with the standard analysis chain described in Chp. 5. 
The extended cuts (Tbl. 5.3) have been applied with both the ring-background and region- 
background model (Sec. 5.3.3). The ON-region was chosen near the centroid of the excess 
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Figure 8.2: Zenith angle distribution of 
7-ray candidates from the H.E.S.S. data 
for MSH 15-52. The peaks indicate the 
different zenith angles of culmination of 
MSH 15—52 for each month of the obser- 
vation period. 
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(15'^14™4?8, -59°9'36") as obtained from a fit of a Gaussian function (Tbl. 8.4). With a radius 

of 0.3°, the ON-region includes most of the y-ray events. To avoid systematic errors from bright 
stars, the OFF-regions have been chosen to exclude stars with a magnitude brighter than seven. 
Also, the ON-region is free of bright stars and therefore this configuration of the background 
models provides good conditions for a reliable analysis. 

8.2.1 Ring-Background Model 

Fig. 8.3 shows the excess map of MSH 15—52 which was obtained with the ring-background 
model. ON- and OFF-regions as applied here are shown by the white circles. The ring radii are 
0.725 and 1.075°. Stars with a magnitude brighter than 7.8 are shown with the labels of their 
magnitudes. The ring-background yields a clear detection of a y-ray excess with a significance 
of 28 standard deviations according to Li and Ma [1983]. The corresponding significance map 
is shown in Fig. 8.4. 

The analysis was repeated with the same configuration but an increased minimum image 
amplitude cut of 400 p.e. This cut implies a higher energy threshold of about 900 GeV and 
therefore a higher reconstruction accuracy and smaller PSF. Although this cut reduces the num- 
ber of excess events to ~25%, the signal to noise ratio and the significance is increased. The 
400p.e. map reveals a slightly different picture of MSH 15—52 (Fig. 8.15), which is discussed 
in Sec. 8.6. The statistics of both analyses are summarized in Tbl. 8.3. 

8.2.2 Region-Background Model 

The analysis was also repeated with the region-background model. The configuration of the 
ON- and OFF-regions is shown in the significance map of Fig. 8.4, which was produced by the 
ring-background described in the previous section. One OFF-region was chosen for each ob- 
servation position opposite to the ON-region, resulting in a total of four different OFF-regions. 
The restriction to one OFF-region was consequence of the partially small distances from the 
ON-region to the four wobble positions. To exclude a star of magnitude 4.5, the OFF-region 
to the left was shifted in a clockwise direction along the arc of constant acceptance. Its previ- 
ous position is shown by the dashed circle. In the final configuration no stars brighter than a 
magnitude of 7.3 are included in the ON- and OFF-regions. Again, a clear detection of a y-ray 
signal is found with a significance of 23 standard deviations. The statistics are given in Tbl. 8.3. 
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Figure 8.3: Excess map obtained with the ring- 
background model showing the ON- and OFF- 
region by white circles. Bright stars are in- 
dicated. They lie outside the ON- and OFF- 
regions. PSR B1509— 58 is represented by the 
black circle. 



Figure 8.4: Significance map obtained with the 
ring-background model. It shows the ON- and 
OFF-regions of the region-background model 
(white circles). Bright stars are indicated. The 
dashed circle shows the position of the left OFF- 
region before its shift to a new position to ex- 
clude a star of magnitude 4.5. PSR B 1509-58 
is represented by the black circle. 



These results are in good agreement with the results of the ring-background model. The same 
configuration was used for the determination of the energy spectrum. 

Table 8.3: Signal statistics as obtained with the standard analysis and the extended cuts (Tbl. 5.3) 
for different background models and image amplitude (M) cuts. 



Ring-background Region-background Ring-background 
(M>80p.e.) (M>80p.e.) (/A>400p.e.) 



A^ON 1V651 17760 1333 

A^OFF 74895 13711 2078 

a 0.186 1 0.196 

Ny 3752 ±133 4049 ±177 925 ±38 

S[o] 27.9 22.9 31.6 

S/Vt[<y/Vh] 5.46 4.47 6.18 

signal /noise 0.27 0.30 2.27 

7-rate [min-i] 2.40±0.09 2.58±0.11 0.59±0.02 
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8.3 Position and Size of the y-ray Excess 



MSH 15—52 is located close to the galactic plane at a distance of 1° as shown in Fig. 8.1. The 
position and the intrinsic extension of the y-ray excess were determined from a fit of the Gaus- 
sian function ,j,,,a(),5o,ffl,A'(^' ^) excess map of Fig. 8.3 as described in Sec. 5.3.5. The 
explicit formula for G is given in Eqn. C.6 of App. C. The PSF fit parameters where chosen ac- 
cording to parameterization number 3 of Tbl. 5.1. The ;f^/dof. of the fit was 391.2/394, equiva- 
lent to a probability of 0.53 for the approximation of the y-ray excess by a Gaussian model. The 
best fit position of the emission was found at (J2000) (15^14™6.'5 ± 2.H, -59°10'1.2" ± 21") 
and the H.E.S.S. catalog name HESS J1514— 591 was assigned to it. The distance to the radio 
position of PSR B1509-58 (15'^13'"55.'620, -59°08'09f'0) is 2.3'. The galactic coordinates of 
the fit position are (320.324° ± 0.005°, -1.200° ± 0.005°). The emission region is asymmet- 
ric and elongated in the north-west direction, defining the major and minor axis with a length 
and width of 6.5' ±0.5' and 2.3' ±0.4' respectively. The angle between the major axis and the 
RA-axis is 43° ±4°. Tbl. 8.4 summarizes the fit results. Fig. 8.5 shows the excess map with 
the contour lines of the fit function (white, solid) and of the intrinsic Gaussian (green, dashed) 
at levels of 25, 50 and 75%. Fig. 8.6 shows slices along the major and minor axes of the excess 
map. The fit function and its components, i.e. intrinsic width or length and the PSF, are shown 
in addition to the data. 
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Figure 8.5: y-ray excess map of the MSH 
15—52 region. Contour lines at levels of 25, 50 
and 75% indicate the fit function (white, solid) 
and the function of the intrinsic width (green, 
dashed). The black cross represents the best fit 
position and the statistical error. The systematic 
error has about the same size. The position of 
PSR B 1509-58 is indicated by the black circle, 
a profile of the PSF is shown in the bottom right. 
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Figure 8.6: y-ray excess counts and fit func- 
tions along the major axis (upper plot) and 
minor axis (lower plot) of the y-ray excess. 
The fit function as well as the functions of the 
intrinsic width and of the PSF are shown. 



105 



8.4 Energy Spectrum 



CHAPTER 8 Detection ofMSH 15-52 



Table 8.4: Parameters obtained from the fit of a Gaussian function (Eqn. C.6) to the y-ray excess 
map of MSH 15—52. The parameterization of the PSF is given in Tbl. 5.1. 



Parameter 


Value 


RA (ao) 


IS'^M'^e.^S ± 2.M, (228.527° ± 0.010°) 


Dec (5o) 


-59° 10' 1 '.'2 ±21", (-59. 169° ± 0.005°) 


Length {a^) 


6.5' ±0.5', (0.109° ±0.008°) 


Width (Oy) 


2.3' ±0.4', (0.039° ±0.006°) 


Angle with RA-axis (©) 


43° ±4° 


PSF parameterization no. 


3 



8.3.1 Systematic Errors 
Position 

The systematic error of the position is determined by the pointing accuracy of the telescope 
system, which is 20" in each direction. In RA the error scales as 1/ cos (Dec) and the relation 
between a second of arc to second of RA is 15"=1 s. So at the declination of MSH 15—52 the 
error in RA is 20"/ cos(— 59. 16°) x 1 s/15" = 2.6 s. Including systematic errors, the fit position 
of the excess is (J2000) (15hl4'"6.^5±2.Hstat±2.'6syst,-59°10'lf'2±21'/tat±20^'y3t). 

Size 

A systematic error of the width and length of the y-ray excess can result from an imperfect 
model of the PSF, which is what was used in the fit function (G, App. C). Fig. 5.6 shows that 
the width of the PSF is slightly overestimated by the Gaussian parameterization, which would 
result in an underestimation of the intrinsic width and length of a y-ray excess. One can reduce 
the PSF parameter 02 by 20% to obtain a model of the PSF with a similar accuracy but slightly 
underestimated width. In this case, the intrinsic width and length of the y-ray excess of MSH 
15—52 would be overestimated by ~5% and ~1% respectively. These numbers provide an 
estimate for the systematic uncertainty in width and length of the y-ray excess. With respect to 
the statistical errors they are negligible. 

8.4 Energy Spectrum 

The energy spectrum was determined using the extended cuts and the region-background model 
with the configuration described in Sec. 8.2. The effective area (Aexi) was obtained from simu- 
lations of an extended source with a standard deviation of 0.04° and 0.1 1° in width and length 
respectively (cf. Sec. 5.4). The spectral parameters were determined by a ;t^-fit to a power law 
according to Eqn. 5.20. It provides a good description for the energy spectrum with a /dof 
of 13.4/12 corresponding to a probability of 0.34. A photon index (F) of about 2.32 ± 0.04 and 
a differential flux at 1 TeV (^ixev) of (5.8 ±0.2) x 10 ^^cm iTeV ^ were obtained. The 
integrated flux above 1 TeV {^{E > ITeV)) is (4.4 ±0.2) x 10^^^cm^^s^\ corresponding to 
about 20% of the flux from the Crab Nebula above the same threshold. The safe energy thresh- 
old was found at ~280 GeV, allowing for a fit range from 280 GeV to 40 TeV. The results are 
summarized in Tbl. 8.5. The data and the fit function are shown in Fig. 8.7. 
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Figure 8.7: Energy spectrum and power law fit of the /-ray emission from the region of MSH 15—52 
with radius 6 = 0.3° in the energy rang from 0.28 to 40 TeV. 



Fig. 8.8 shows the one, two and three (7 error contours of ^ixeV and T as obtained from the 
fit with the TMinuit class in ROOT. The contours correspond to a ;f^/dof of 2.3/2, 6.18/2 and 
11.83/2 with the probabilities of 0.68, 0.95 and 0.99 respectively. 

To investigate the emission region for spatial variations of the photon index, the analysis 
was repeated for two different smaller signal regions with radius 6 = 0.045° — one near the 
center of the emission at the position of PSR B 1509— 58 and the other one in some distance at 
the X-ray jet axis to the southeast. The regions are shown in Fig. 8.9. The area and position 
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*{cm'^s-'TeV') 

Figure 8.8: One, two and three a error contours 
for the spectral parameters F and 0iTeV as ob- 
tained from the ;^^-fit to a power law in Fig. 8.7. 
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Figure 8.9: Excess map showing the y-ray ex- 
cess region (6 = 0.3°, white) and two smaller 
ON-regions (red and green) where the energy 
spectrum was determined. The position of PSR 
B1509— 58 is indicated by the black circle. 
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of the background regions, and the effective area, were adjusted accordingly. The measured 
energy spectrum was fit to a power law in the same fit range, but no significant variation of the 
photon index was found. In the pulsar region, a F of 2.22 ± 0.08 with a ;^^/dof of 5.0/6 and in 
the jet region, a F of 2. 16 ± 0. 10 with a ;^^/dof of 7.0/5 was obtained. For the interpretation 
of the results one has to consider that the regions are slightly correlated due to the PSF and 
that the ON-regions are slightly offset from the center of the excess. While the first effect 
reduces spectral differences between the regions, the later leads to a spectral steepening. This is 
explained by the energy dependence of the PSF, which spreads low energy events further from 
the center. However, the effect on the spectral index is very small. 

8.4.1 Systematic Errors 

Since the energy spectrum is determined at a late point in the analysis chain, it accumulates 
errors from previous analysis steps. Therefore it is important to consider some systematic un- 
certainties. 

Atmospheric Model of the Simulations 

The calculation of the energy spectrum is sensitive to atmospheric variations. Unfortunately 
it is very difficult to model the atmospheric variations in all detail. Therefore, Monte Carlo 
simulations only contain a model of the mean atmospheric density profile. The desert and 
the maritime atmospheric models (Sec. 5.1) are both appropriate models of the atmosphere 
at the H.E.S.S. site, as confirmed with trigger studies by Funk et al. [2004]. In the analysis 
the effective area for the desert model was used. The effective area for the maritime model 
provides slightly different results and these are summarized in Tbl. 8.5. The differences provide 
an estimate for the systematic error as resulting from uncertainties in the atmospheric model. 
They are 15% and 3.5% for 0iTev and F respectively. 

Absolute Calibration 

The degradation of the absolute photon efficiency of H.E.S.S. and its measurement with muons 

was previously discussed in Sec. 5.2.1. Fig. 8.10 shows the degradation of the muon efficiency 
during the observation of MSH 15—52 as determined by Bolz [2004b]. It can be seen that the 
mean muon efficiency during this period is only 87% of the nominal value of 0.106 and also a 
decrease by 5% from the beginning to the end of this period is observed. Further detailed studies 
by Bruno Khelifi and Conor Masterson [2005] for this observation period have shown that in the 
standard analysis the differential flux d^/dE is underestimated by about 13%, while the photon 
index is not affected significantly (Bruno Khelifi and Conor Masterson [2005]). These results 
provide an estimate of the systematic error of the spectrum related to the absolute calibration. 

Cut Configuration 

The cut configuration also contributes a small systematic uncertainty to the energy spectrum. 
To estimate this uncertainty, the analysis was repeated with loose and hard cuts (cf. Tbl. 5.3). 
The results in Tbl. 8.5 show that (^ixeV is either slightly over or underestimated. The systematic 
error of 0iTeV and F is estimated as 7% and 1%. 
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Efficiency vs Timestamp 




Figure 8.10: Muon efficiency during the obser- 
vation period of MSH 15—52. The time is mea- 
sured in units of years since 2000. The aver- 
age efficiency is only 87% of the nominal value 
0.106. The decrease in the efficiency is about 
5% during the observation period. (Figure taken 
from Bolz [2004b].) 
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Broken Pixels 

The broken pixels discussed in Sec. 5.2. 1 introduce a systematic uncertainty of about 5% to each 
image amplitude, as detailed studies by Aharonian et al. [2006a] and Schwanke et al. [2004] 
have shown. This value can be directly adopted as the systematic error of the flux normalization 

01TeV- 

Uncertainty of Sources Extension 

The effective areas were derived from Monte Carlo simulations for extended sources as de- 
scribed in Sec. 5.1.3. The standard error of the extension is about 0.07° (cf. Tbl. 8.4) in 
width and length. To determine the resulting uncertainty for the energy spectrum, the spectrum 
was reproduced with effective areas for an increased (o^, — 0.046°, (7/ = 0.12°) and reduced 
(Ow = 0.034°, Oi = 0.10°) source extension by approximately one standard error. Tbl. 8.5 
shows the results of a power law fit to the two new energy spectra. The relative systematic 
uncertainty for ^ixeV is about 2%, while T is unaffected. 

Simulation of Extended Sources 

A marginal systematic error results from the simulation of extended sources which do not take 
the radial acceptance gradient in the FOV into account. Fig. 8.11 shows the radial profile of the 
system acceptance. The acceptance decreases by about 1% per 0.1°. The estimated total event 
loss for a source with a size of MSH 15—52 is roughly 1%, providing also an estimate for the 
error of the ^ixeV- 
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1 
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Live-time 

There is an uncertainty of 1% in the determination of live-time. This was determined by Funk 
et al. [2004]. The error of the live-time translates directly into an error of the flux normalization. 



Total Systematic Error 

Tbl. 8.5 shows the spectra as obtained for the different analysis configurations discussed above. 
Tbl. 8.6 shows the relative systematic errors of O and P. Under the reasonable assumption that 
the individual systematic errors are uncorrelated, the total relative systematic error is obtained 
from the quadratic mean of individual systematic errors. One finds A0iTev/0iTeV = 22% and 
Ar/r = 4%. The total relative systematic error of the integral flux (A^/^))(£' > ITeV) = 23% 
is obtained by error propagation according to Eqn. 5.23. The covariance is zero in this case, 
since the errors in 4> and T are assumed to be independent. 



Table 8.5: Energy spectra as determined from a fit to a power law for the different configurations 
discussed in this section. The meanings of superscript numbers are as follows: *' configuration 
providing the quoted results, ^with maritime atmospheric model, with hard cuts, ^ with loose 
cuts, with decreased (csw = 0.034°, CJW = 0.10°) source extension and ^ with increased (o^ = 
0.046°, = 0.12°) source extension. 



Configuration 


01TeV 


r 


^{E >lTeV) 


;t:^/dof 






[lO-'^cm-^s-'TeV-i] 




[lO-i^cm-^s-i] 






extended " 


5.84±0.23 


2.32±0.04 


4.43±0.18 


13.4/12 


0.34 


maritime atm ^ 


6.74±0.26 


2.24±0.04 


5.46±0.22 


10.4/8 


0.24 


hard cuts 


6.20±0.27 


2.34±0.04 


4.62±0.17 


22.7/20 


0.30 


loose cuts ^ 


5.36±0.28 


2.29±0.04 


4.14±0.20 


10.6/11 


0.48 


4 

size- 


5.76±0.23 


2.32±0.04 


4.37±0.17 


13.4/12 


0.34 


size+ 


5.97±().24 


2.32±().()4 


4.52±().18 


13.4/12 


0.34 



Table 8.6: Individual and total relative systematic errors for the energy spectrum. 



Systematic Error A0iTev/<^iTev[%] Ar/r[% 



Atmospheric model 15 3.5 

Absolute calibration 13 

Cut configuration 7 1 

Broken pixels 5 

Uncertainty of source extension 2 

Simulation for extended sources 1 

Live-time 1 

Total 22 4~ 
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8.5 Light Curve 

A light curve of the daily integrated flux average ^»(£' > ITeV) for the emission from the ON- 
region (0 = 0.3°) is shown in Fig. 8.12. It is obtained from a power law fit to the daily data 
assuming a constant photon index of F = 2.3. A fit of the light curve to a constant (red line) 
yields a 4>(£ > ITeV) of (4.1 ±0.2) x lO^^^cm-^s"! and a ;i:^/dof of 26.2/21. This resuk is in 
good agreement with a constant y-ray flux from MSH 15—52. 



Figure 8.12: Light curve showing the 
daily integrated flux average <!>(£'> 1 TeV) 
from the region of MSH 15-52. The fit to 
a constant (red line) is in good agreement 
with a constant 7-ray emission. 
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8.6 7-Ray Morphology 

The morphology of MSH 15—52 is difficult to resolve with H.E.S.S., since the extension of the 
7-ray excess is small and of the same order as the size of the PSF Nevertheless, a better knowl- 
edge of the 7-ray excess distribution is desirable for a better picture of the processes of the PWN 
and for a comparison with structures resolved in X-rays. Here an attempt is made to investigate 
the 7-ray morphology of MSH 15—52 in greater detail using different representations, energy 
bands and methods. 



8.6.1 Sky Maps 

Fig. 8.13 to 8.16 shows four different 7-ray maps which differ in their image amplitude cuts and 
methods of restoration. An increased image amplitude cut provides an increased resolution at 
the expense of an increased energy threshold. Deconvolution (cf. Chp. 6) on the other hand can 
also reduce the apparent PSF while preserving the original energy range. However, statistical 
artefacts are amplified. The method of image smoothing can always provide a conservative 
representation. The sky maps are based on the count map described before. A profile of the 
PSF is shown in the bottom right of each map. 

Fig. 8.13 shows the smoothed excess map with a minimum image amplitude cut of 80p.e. 
with an energy threshold of about 280 GeV. The smoothing Gaussian has a width Os = 0.03°. 
An elliptical 7-ray excess region with small deviations from a Gaussian distribution can be seen. 
PSR B 1509— 58 (black circle) lies close to the peak intensity. The containment radius of 68% 
is 0.12° (Sec. 5.2.6). 

Fig. 8.14 shows the deconvolved count map. It was obtained using the same cuts and 10 
iterations with the Richardson-Lucy algorithm. The results are similar to the smoothed excess 
map of Fig. 8.13, but the overall width is smaller and the peak intensity is higher since the 
spread by the PSF is reduced. To make Fig. 8.14 comparable to Fig. 8.13, the background of 
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Figure 8.13: y-ray excess map convolved with 
a Gaussian (a = 0.03°). The position of PSR 
B 1509— 58 is marked by the black circle. The 
PSF is indicated in the bottom right corner. 
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Figure 8.14: 7-ray count map after 10 iterations 
with the Richardson-Lucy algorithm. The posi- 
tion of PSR B 1509-58 is marked by the black 
circle. The PSF is indicated in the bottom right. 
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Smoothed Excess Map (IA>400 p.e.) 
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Figure 8.15: The same as Fig. 8.13, but with a 
minimum image amplitude of at least 400 p.e., 
reducing the width of the PSF by about 50%. 
The PSF is indicated in the bottom right corner. 
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Figure 8.16: The same as Fig. 8.14, but with a 
minimum image amplitude of at least 400 p.e. 
reducing the width of the PSF by about 50%. 
The PSF is indicated in the bottom right comer. 



5.0 counts/bin has to be subtracted. Then the peak intensities of the smoothed and the restored 
map compare as 5 to 9 counts. 

A slightly different image is obtained with an image amplitude cut of 400 p.e., which raises 
the energy threshold to about 900 GeV. Fig. 8.15 shows this excess map smoothed with a Gaus- 
sian function (c7 = 0.03°) . Due to the increased energy threshold, the size of the PSF is reduced 
to about half its size at 80 p.e. increasing the resolution. The 68% containment radius is 0.063° 
(Sec. 5.2.6). Although the event statistic is reduced, the signal to noise ratio is higher. The main 
emission region in this map appears slightly curved. 
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Fig. 8.16 shows the 400 p.e. count map after deconvolution with 10 iterations of the Richard- 
son-Lucy algorithm. The emission region has a significantly reduced width but about the same 
length as the corresponding smoothed map. The highest intensity is seen in a compact region 
oriented along a northwest- southeast direction which covers the X-ray jet axis. The emission 
region is also slightly curved and extends to the southwest of the pulsar. 

The maps are interesting and complementary to each other: the 80 p.e. maps for aspects 
to a low energy threshold, the 400 p.e. maps for their increased resolution. The also suggest a 
slightly different morphology at low and high energies motivating an investigation for an energy 
dependent morphology in the next section. 

8.6.2 Energy Bands 

To investigate the energy dependence of the morphology, the data was divided in five differ- 
ent energy bands with the energy ranges 0.2-0.5, 0.5-1, 1-2, 2-5 and 5-100 TeV. The event 
statistics and PSF of each energy band is summarized in Tbl. 8.7 and Tbl. G.l respectively. 
The corresponding excess maps are shown in Fig. 8.17 (left column) each with a profile of the 
PSF in the bottom right. The intrinsic extension of the 7-ray excess of each energy band was 
determined from a fit of the Gaussian function (G) to the map as described in Sec. 8.3. The 
contour lines of G (white) and the component representing the intrinsic source size (green) are 
overlaid. One-dimensional slices of the fit functions are shown in App. G. 

Fig. 8.18 show the intrinsic standard deviations of the major and minor axis versus the 
different energy bands in order of increasing energy (cf. Tbl. 8.7). A decreasing trend of the 
longitudinal source extension is found. Fig. 8.19 and 8.20 show the same representation for 
the centroid of the fit function in RA and Dec. In both cases the centroid converges towards 
the pulsar position with increasing energy. The pulsar position is indicated by the horizontal 
dashed lines. This is what would be expected if the pulsar is the source of a y-ray producing 
pulsar wind. 

The second column in Fig. 8.17 shows the excess maps smoothed with a Gaussian function 
of standard deviation Os = 0.03°. Again the size of the smoothed emission region decreases 
with energy, but a immediate conclusion about the intrinsic source extension is complicated by 
the energy dependence of the PSFs (cf. Sec. 5.2.6). 

However, a conclusion becomes possible, if the PSF is similar for each energy band. There- 
fore count maps with a smaller PSF, have been artificially smoothed to increase the PSF to the 
same 68% containment radius of 0.125°. The difference of the source extensions between high 
and low energies is more apparent if the sky maps below 5 TeV are combined. The required 
smoothing values {<Js) are 0.023° and 0.064° respectively (cf. Tbl. G.l). The resulting two sky 
maps are shown in the third column of Fig. 8. 17. 

Fig. 8.21 shows the combination of these two smoothed maps. Two complementary color 
scales, i.e. yellow for E <5 TeV and blue for £ > 5 TeV, have been chosen to provide a color 
neutral, i.e. white, appearance for similar bright regions. However, the color brightness between 
both energy bands is not normalized. Therefore, it does not represent absolute numbers of events 
but the relative intensity profile of each band. The 68% and 95% containment radii are shown 
at the bottom left. The later do not match exactly since the PSF parameterization changes with 
energy and therefore prevents this. Since the central region appears white and the outer region 
yellow, the 7-ray emission region is more compact at higher energies. 
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Excess Map, 0.2TeV<E<0.5TeV 



Smoothed Excess Map, 0.2TeV<E<0.5TeV 
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Smoothed Excess Map, 0.5TeV<E<1TeV 
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Smoothed Excess Map, 1TeV<E<2TeV 
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Smoothed Excess Map, 2TeV<E<5TeV 



Figure 8.17: Morphology at 
different energy bands. 
First column: excess map 
overlaid with Gaussian fit 
function (white) and the cor- 
responding component of the 
intrinsic width (green). 
Second column: smoothed 
excess map ((7s=0.03°). 

Third column: excess maps 
additionally smoothed to the 
same 68% containment radius 
of 0.125°. The sky maps for 
E < 5TeV have been com- 
bined. The PSF is indicated 
in the bottom right of each 
plot. The position of PSR 
B 1509-58 is marked by the 
black circle. 
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Excess Map, 5TeV<E<100TeV 



Smoothed Excess Map, 5TeV<E<100TeV 
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Figure 8.18: Intrinsic width 
and length of the Gaussian 
function G (App. C) for the fit 
to the excess maps of the en- 
ergy bands in Fig. 8.17. 
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Figure 8.19: Centroid in RA 
of the Gaussian function G 
(App. C) for the fit to the ex- 
cess maps of the energy bands 
in Fig. 8.17. 
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Figure 8.20: Centroid in Dec 
of the Gaussian function G 
(App. C) for the fit to the ex- 
cess maps of the energy bands 
in Fig. 8.17. 




Figure 8.21: Two color image of 
MSH 15—52 combining the maps from 
Fig. 8.17 for E < 5TeV (yellow) and 
E > 5TeV (blue). The color palettes 
have been chosen complementary to each 
other to provide a color neutral (white) 
appearance for similar brightness. How- 
ever, the color brightness between both 
energy bands is not normalized. There- 
fore it does not represent the absolute 
number of events but the relative inten- 
sity profile of each band. The position of 
PSR B 1509-58 is indicated by the circle. 
The 68% and 95% containment radii are 



15h16m i5h14m 15h12m shown at the bottom left, the profile of the 

RA (hours) ^ 

PSF at the bottom right. 



Table 8.7: Statistics and significance {S) for different energy bands as obtained from the ring- 
background model with the same configuration discussed in Sec. 8.2.1 before. The energy bands 
have been numbered in oder of increasing energy. 



No. 


E [TeV] 


^ON 


A^OFF 


Ny 


a 


S[a] 


1 


0.28 TeV<£< 0.5 TeV 


10820 


51594 


mi 


0.185 


11.8 


2 


0.5TeV<£ < ITeV 


3615 


13165 


1051 


0.195 


17.7 


3 


1 TeV< E <2 TeV 


1875 


5532 


814 


0.192 


20.3 


4 


2 TeV< £ < 5 TeV 


1045 


3932 


422 


0.158 


14.1 


5 


5TeV<£ < 100 TeV 


296 


672 


185 


0.165 


13.0 




0.28TeV<£ < 100 TeV 


17651 


74895 


3752 


0.186 


27.9 



In summary, the energy dependent analysis supports a picture of a decreasing longitudinal 
extension of the y-ray emission along the pulsar jet axis with increasing energy. This is similar 
to what is observed at X-rays (Forot et al. [2006]) as pointed out in Sec. 2.1.1. 
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8.6.3 Correlation with X-ray Emission 

The similarity between the y-ray and X-ray morphology discussed in Chp. 2 motivates further 
comparisons. According to Chp. 2.4.5 a correlation between X- and y-ray emission is expected 
for a synchrotron and IC radiation producing wind of VHE electrons. Here this correlation is 
investigated using H.E.S.S. 7- and X-ray data from the ROSAT and Chandra satellites. 

ROSAT Data 

The ROSAT X-ray data was recorded in 1991 and 1992 covering the energy range from 0.1 
to 2.4keV. Fig. 8.22 shows the deconvolved H.E.S.S. 7-ray map of Fig. 8.16 overlaid with the 
contour lines of the ROSAT data at the levels of 0.5, 1, 1.5, 2, 5, 10, 20, 30 and 40 in arbitrary 
units as determined by Trussoni et al. [1996]. A good correlation between the 7- and X-ray 
data is found. However, it is interesting to note that the region of RCW 89 shows strong X-ray 
emission but no strong 7-ray emission. Such a difference can indicate a difference of the X-ray 
production mechanisms in this region. 



Figure 8.22: ROSAT X-ray 
contours overlaid on the recon- 
structed H.E.S.S. 7-ray map of 
Fig. 8.15. A good correlation 
between both wavelengths can 
be seen. A profile of the H.E.S.S 
PSF is indicated in the bottom 
right comer. The position of 
PSR B 1509-58 is marked by 
the black circle. 




15h16m 15h14m 15h12m 
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Chandra Data 

In 2000, the Chandra satellite provided new X-ray data of higher resolution and statistics. The 
data discussed here was recorded on August 14, 2000 with observation ID 754. It is the same 
data which was discussed in Chp. 2 based on the analysis carried out by Gaensler et al. [2002]. 
It has a live-time of 19039 s, a total number of 222914 X-ray events and covers the energy 
range from 0.3 to lOkeV. Fig. 8.23 shows the exposure-corrected X-ray map of this data on a 
logarithmic scale. Details about the Chandra data analysis and exposure correction are given 
in App. H. The Chandra and the ROAST X-ray maps are in good agreement. Both show 
significant X-ray emission from the region of PSR B 1509— 58 and RCW 89 and an elongated 
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Chandra Exposure Corrected Map 
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Figure 8.23: Chandra map of the X-ray emis- 
sion in the energy range from 0.3 to lOkeV. 
Strong X-ray emission is seen from the region 
of PSR B 1509-58 (black circle) and the north- 
west region of RCW 89. 
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Figure 8.24: H.E.S.S. 7-ray excess map 
overlaid with the Chandra X-ray contours of 
Fig. 8.23. The contour lines are chosen at the 
levels of 3, 6, 12, 18, and 30% of the peak inten- 
sity. The PSF is shown (bottom right). 



structure extending to the southeast. Moreover, the Chandra map can clearly resolve the jet- 
like character of this structure and the clumps of X-ray emission in the region of RCW 89. 
Fig. 8.24 shows the contour lines of the Chandra data at levels of 3, 6, 12, 18 and 30% overlaid 
on the H.E.S.S. image of Fig. 8.13. The contour lines were obtained from the count map after 
smoothing with a Gaussian function of a width of 0.005°. The Chandra field of view is indicated 
by the white square. Again, the y-ray emission seems to be correlated with the X-ray emission 
at the PSR B 1509-58 region, however not at the RCW 89 region. 

To quantify the correlation between the 7- and X-ray emission, the correlation coefficient 
was calculated for the H.E.S.S. and Chandra data in two rectangular regions. One of them 
(^pulsar) covers the region of PSR B 1509— 58 and the other (i?Rcw89) covers the region of 
RCW 89. The correlation coefficient (p) is defined as 

, , ^l'ir.--rm-x) ^^^^ 
v'If.i()i-r)-If.i«-^fP 

where 7 and represent the content of bin / of the 7- and X-ray maps respectively. 7 and 
X denote the corresponding mean values, p is bound to the interval (—1 < p < 1), where 
the values 1, 0, and -1 express a perfect correlation, no correlation and perfect anti-correlation 
respectively. 

The errors of the correlation coefficient (Op, Op) due to errors in 7 (Oyj) and Xt {ox,i) 
are asymmetric, as implied by the asymmetric probability density distribution of p on the fi- 
nite interval. Therefore they cannot be immediately calculated by Gaussian error propagation. 
However, Fisher's Z- Transform (Fisher [1925]) provides the appropriate transformation of p to 
the Gaussian probability density distribution ^(p) on the infinite interval (—00 < CiP) < 
where Gaussian error propagation can be applied and Gaussian confidence intervals are valid. 
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Fisher's Z- Transform (^) is defined as 

C(p) = atanh(p) (8.2) 

and has the following properties 

C(p-ap-) = C(p)-(T^(p), (8.3) 

C(p + (Jp+)-C(p) + af(p), (8.4) 

where is the symmetric error of ^. Hence with the inverse transformation (^~^ = tanh) the 
asymmetric errors are 

a- = - [i;-' (C(p) - C7^(p)) -p] = - [tanh(atanh(p) -a^(p)) -p] , (8.5) 

a; = + [C"' (C(P) + -P] = + [tanh(atanh(p) + C7^(p)) -p] . (8.6) 

(7^(p) can be found by Gaussian error propagation as 



and 




(8.7) 



(8.8) 



The calculation of the partial derivatives dp/dji and dp/dXi is described in Schwanke and 
Lohse [2005]. For the H.E.S.S. data C7y,,- = a/A^onj + «^A'^OFF.i- For the Chandra data Oxd = 
\/Xi, where X is the number of X-ray counts before the exposure correction. Fig. 8.25 shows 
the y-ray excess and X-ray count map with equal binning as well as with the two regions where 
the correlation coefficient has been determined. The bin size of the maps has been increased to 
0.02° X 0.02°. The /-ray count map has been convolved with a Gaussian of 0.02° before the 
equally convolved background map was subtracted. The regions i?puisar and i?Rcw89 consist of 
9x8 and 6x4 bins in RA and Dec respectively. The coordinates and which define 
the lower left comer of each region, are given in Tbl. 8.8. The corresponding scatter plots of the 
7" and X-ray maps are shown in Fig. 8.26. The counts in each bin are sufficiently high enough 
to justify the assumption of Gaussian errors. 

At the regions of i?puisar and i?Rcw89 correlation coefficients of 0.58^q|^ and — 0.06^q33 are 
obtained respectively. They support the speculations about a difference of the X-ray production 
mechanisms in the region of RCW 89 in comparison to the region of PSR B 1509— 58. 
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15h15m 15h14m 15h13m 15h15m 15h14m 15h13m 
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Figure 8.25: H.E.S.S. /-ray excess map (left) and Chandra X-ray map (right) of MSH 15—52 with 
identical binning and a bin width of 0.02° x 0.02°. The two adjacent (white) squares define the 
regions /?puisar (9 bins x 8 bins) and Rrcw&9 (6 bins x 4 bins) where correlation coefficients are deter- 
mined. The position of PSR B1509— 58 is indicated by the black circle. 




Y-ray counts Y-ray counts 



Figure 8.26: Scatter plots for determining the correlation coefficient between the H.E.S.S. 7-ray and 
the Chandra X-ray data in the regions of PSR B 1509-58 ( /?puisar) and RCW 89 (/?i?cw89) which are 
indicated in Fig. 8.25. The peak value is outside the visible range of the X-ray axis in each of the 
plots. 



Table 8.8: Correlation between the H.E.S.S. 7-ray and Chandra X-ray data in the regions of PSR 
B 1509— 58 (/?puisar) and RCW 89 {Rrcw^9)- The coordinates mark the lower left corners of each of 
the two square regions shown in Fig. 8.25. 





^pulsar 


'^RCW89 




15'^14'"35.'23 


15'^14™7.n4 




-59.245° 


-59.085° 


No. of bins (RAxDec) 


9x8 


6x4 


P 


0.58l«:|3 


-0 06+^-34 

^•^"-0.33 
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8.7 Search for Pulsed Emission from PSR B1509-58 

Since a strong y-ray signal from the region of PSR B 1509— 58 was detected it is possible and 
interesting to analyze this signal for pulsed emission with a period of the pulsar. In contrast 
to the constant emission seen from the extended PWN of MSH 15—52, the detection of pulsed 
emission would indicate emission processes near the light cylinder of PSR B 1509— 58. Since 
this has never been observed for any pulsar at the VHE range accessible to H.E.S.S. it would 
have been a very remarkable finding. In order to test for pulsed emission, the data was analyzed 
for periodicity according to the radio ephemeris of Tbl. 8.9, which was taken from the ATNF 
Data Archive (Manchester et al. [2006]). The extended cut configuration was chosen for this 
analysis and applied to a reduced ON-region shown in Fig. 8.9 (red circle). The small radius 
(6 = 0.045°) reduces the number of y-rays events which do not originate near the pulsar. The 
events in the ON-region have been filled into a phasogram as described in Chp. 7. The resulting 
phasogram in Fig. 8.27 shows a uniform distribution. It covers the full H.E.S.S. energy range 
from 280 GeV to 40TeV. 

Since most pulsar models predict a decrease of the pulsed VHE flux with increasing energy, 
the detection of pulsed emission is more likely to occur at the lower end of the energy spectrum. 
Consequently, removing events of higher energy reduces events which are less likely pulsed and 
therefore increases the ratio of pulsed y-ray flux — at a cost of event statistics. So the analysis 
was repeated for a reduced energy range of 280 to 500 GeV. Fig. 8.28 shows the corresponding 
phasogram. Also this light curve is consistent with a uniform distribution. 

Table 8.9: Parameters of the ephemeris of PSR B 1509-58 from the ATNF archive (Manchester 
et al. [2006]) which have been used in the periodic analysis. 



Parameter Value 
RA (J2000) 15^13'"55.^620 
Dec (J2000) -59°08'9f'00 
validity range [MJD] [53035; 53554] 

t^'^^ [MJD] 53295.000000712 

/o [s-i] 6.6088537688620 

/o[s-2] -6.68672x 10 11 



8.7.1 Tests for Periodicity 

The uniformity of the distribution and the non-detection of pulsed emission was confirmed using 
the Pearson test, the Z test and the H test as described in Chp. 7. The results for the reduced 
and full energy ranges are summarized in Tbl. 8.10. Within statistical errors, no indication for 
a pulsed y-ray emission is found. 
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Phasogram of PSR B1509 (280GeV<E<40TeV) 




Phasogram of PSR B1509 (280GeV<E<500GeV) 



Figure 8.27: Phasogram of the PSR B 1509-58 
in the energy range from 0.28 to 40TeV. The 
hatched phase region [0.2; 0.6] was used for the 
calculation of an upper limit of the pulsed flux. 
The mean is indicated by the horizontal line. 




Figure 8.28: Phasogram as described in 
Fig. 8.27 but for the reduced energy range of 
0.28TeV<£ <0.50TeV. 



Table 8.10: Results for the x^, -Z^ and H test (cf. Chp. 7) as applied to the phasograms of Fig. 8.27 
and 8.28. The corresponding probabilities (P) are also listed. They are in agreement with a uniform 
distribution, i.e. no pulsed signal is identified. 



E [TeV] 


x' 




^2 




H 


P{H) 


0.28 <E<40 


24.7 / 24 


0.42 


4.18 


0.38 


3.94 


0.21 


0.28 <£< 0.5 


24.4 / 24 


0.44 


7.68 


0.10 


7.51 


0.05 



8.7.2 Flux Upper Limit 

An upper limit for the pulsed integrated flux was calculated in the phase region selected by 
Kuiper et al. [1999], who reported the detection of pulsed emission from PSR B1509— 58 at 
the energy range of 0.75-lOlVIeV in the phase region [0.2; 0.6] (cf. Sec. 2.1.2, Fig. 2.9). The 
corresponding On- / Off-statistics are given in Tbl. 8.11. The upper limit was determined fol- 
lowing the unified approach by Feldman and Cousins [1998] as described in Chp. 7. The upper 
limit for the integrated flux above the H.E.S.S. energy threshold of 280 GeV and above 1 TeV, 
is shown in the same table for a confidence level (CL) of 99%. The upper limits are given with 
and without the systematic flux error of 22%. 

Table 8.11: Event statistics and the upper limit in the phase region [0.2; 0.6] for the phasogram of 
Fig. 8.27. The upper limits are calculated with (UL^^^) and without {UL(S>) the systematic error for 
a confidence level of 99%. 



Eo [TeV] 


On 


Off 


a 


Excess 


ULct>{E > Eo) 
[10-i2cm-2s-i] 


UL%'\E > Eo) 
[10-i2cm-2s-i] 


0.28 
1.00 


259 
69 


374 
106 


2/3 
2/3 


9.7 
-1.6 


< 8.3 (99% CL) 

< 0.82 (99% CL) 


< 11.0 (99% CL) 

< 0.97 (99% CL) 
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Chapter 9 
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Analysis of the H.E.S.S. TeV y-ray data of MSH 15—52 has provided new infomiation about 
this exotic object. With this information, further conclusions can be drawn when the results are 
interpreted in the context of astrophysical models and together with data from observations at 
other wavelengths. Three such examples are discussed here. They include conclusions about 
the y-ray production in MSH 15—52, about the dominating transport mechanism in the pulsar 
wind and the interaction of the pulsar wind with the optical nebula RCW 89. 

9.1 7-Ray Production in MSH 15-52 

An important discovery from the H.E.S.S. data is that it agrees well with the leptonic model 
of y-ray production in PWN. For example, Aharonian et al. [2005b] showed with simulations 
of GALPROP (Strong et al. [2000]) that the TeV y-ray spectrum can be explained by inverse 
Compton (IC) scattering of cosmic microwave background (CMB) radiation and infrared (IR) 
photons from dust and starlight by highly accelerated electrons. In these simulations, the parent 
energy spectrum of the electrons obeys a power law and is adjusted to reproduce the IC y radi- 
ation measured by H.E.S.S as well as the synchrotron emission determined by X-ray and radio 
measurements by BeppoSAX (Mineo et al. [2001]) and ATCA (Gaensler et al. [2002]). The en- 
ergy density of the ER photons and the magnetic field inside the nebula were kept free to match 
the experimental data. The simulations determine the magnetic field to ~17 [iG, the IR energy 
density to 2.5 eVcm~^ for the IR seed photons with wavelengths of ~100 jim and the spectral 
index of the electrons to Tg ~ 2.9. Fig. 9. 1 shows the corresponding spectral energy distribution 
of the synchrotron and IC radiation, together with the experimental data. The IC radiation from 
the individual components and the sum of the seed photon field of these components is shown. 
Apparently this model provides a good description of the data. 

An alternative interpretation with a model of TeV y radiation from hadronic primary parti- 
cles was introduced in Sec. 2.4.5. Bednarek and Bartosik [2003] have calculated the expected 
y-ray spectrum from MSH 15—52 for nucleonic collisions with a high and low density medium 
surrounding the PWN. The green thick and thin dot-dashed curve in Fig. 9.1 represents the spec- 
tral energy distribution in high and low density mediums with average densities of 300 cm~^ 
and 0.3 cm^-^, respectively. Both curves are taken from Fig. 2.24. It can be seen that in a 
high density medium the predicted TeV y radiation exceeds the observed flux. Therefore the 
H.E.S.S. data constrains the density of the ambient medium below 300 cm~^. On the other 
hand, for the low density medium the expected flux is too low to be visible in the measured TeV 
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spectrum. H I measurements by Dubner et al. [2002] have suggested a density of ~0.4 cm^-^ for 
the interstellar medium in the southeast and a density of ~12cm^ for the northern region, in 
particular a density of ~15 cm^^ for the region of RCW 89. Therefore, in the denser regions, 
in particular the region of RCW 89, one would expect an increased y-ray flux from hadronic 
interactions as pointed out by Bednarek and Bartosik [2003]. The analysis of the H.E.S.S. data 
reveals, however, a decreased y-ray flux in these regions (cf. Sec 8.6). This could imply that the 
nucleon density in the pulsar wind is too low to contribute significantly to the y emission. 

It can be concluded, therefore, that the leptonic radiation model can well explain the ob- 
served TeV spectrum from MSH 15—52. The existence of a significant electron component is 
evident by the synchrotron flux observed in the X-ray band. However, evidence for a hadronic 
y-ray production has not been found and therefore does not play a major role in the TeV y-ray 
production in MSH 15-52. 




10 " 10" 10' 
Energy (eV) 

Figure 9.1: Spectral energy distribution of radiation from MSH 15—52. The radio, X-ray and 
TeV 7-ray measurements by ATCA (Gaensler et al. [2002]), BeppoSAX (Mineo et al. [2001]) and 
H.E.S.S. (Khelifi, B. et al, (H.E.S.S. collaboration) [2005]) are shown along with the simulated 
curves for leptonic and hadronic radiation. The inverse Compton curves are obtained for different 
photon fields. The green thick and thin curves represent the 7 radiation from tt" decay of nucleonic 
interactions in high (300 cm^^) and low (0.3 cm^^) density mediums. (Figure taken from Khelifi, B. 
et al, (H.E.S.S. collaboration) [2005]. Green curves added from Bednarek and Bartosik [2003].) 



9.2 Transport Mechanism in MSH 15-52 

A further, interesting conclusion can be drawn concerning the transport mechanism in MSH 
15—52 if the model by Kennel and Coroniti [1984] is considered. This model assumes that 
particles are highly accelerated near a pulsar and that they drive the expansion of a PWN into 
the ambient interstellar medium (Sec. 2.3.2). Moreover, it assumes that pulsar wind consists of 
a plasma which is frozen into the magnetic field lines of a pulsar. The field lines are wrapped in 
tight spirals around the pulsar, due to the pulsar's spin. Using this picture of pulsar wind, one can 
explain the expansion of this magnetized flow either by diffusion or by convection. Diffusion is 
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the random propagation of scattering particles with a net transport determined by the diffusion 
gradient. Convection, on the other hand, is a collective flow of particles in a common direction. 
According to de Jager [2006c], de Jager [2006a] and de Jager [2006b], one can compare the 
efficiency of these two transport mechanisms in a PWN and determine the dominant transport 
mechanism of the radial flow. 

In the model of Kennel and Coroniti [1984], the pulsar wind is shocked at the shock radius 
rg inside the PWN of radius Rj\f. This situation is illustrated in Fig. 2.15. The flow velocities 
v{Rn) and v{rs) of the wind have to meet the boundary conditions at the radii and Rn, i.e. 
v{Ri\f) is equal to the outer expansion velocity of the nebula and v(rv) = c/^s — c/\/3, where jS^ 
is the expected Alphen velocity of the relativistic plasma (Kennel and Coroniti [1984]). 

Therefore de Jager argued, that a radial flow velocity between and R^ given by 

vW = ^^ (9.1) 

will meet both boundary conditions, and that the 1 /r dependence is a reasonable assumption in 
the case of constant B -field strength (which is generally fulfilled in good approximation) and 
due to conservation of magnetic flux. 

Following the argumentation, one can estimate the diffusive and convective flow in a PWN. 
If the diffusion time is smaller than the convection time Tc, most energetic particles over- 
take the bulk convective flow due to scattering and the transport mechanism is dominated by 
diffusion (and vice versa). Hence the criterion for a flow being dominated by diffusion can be 
written as 

Tc > Td. (9.2) 
The time for diffusion of an electron to a radius r is given by 

W = 5|, (9.3) 

where is the diffusion coefficient for the electron transport perpendicular to the magnetic 
field, i.e. in the radial direction. The diffusion coefficient is given by 

= (9.4) 

where Aj^ is the mean free path between scatterings in the radial direction and v is the velocity 
of the electron. Moreover, in the case of strong and weak scattering the fundamental limit by 
Steenberg [1998] states that 

A±<PL, (9.5) 

where pi is the particle's gyroradius. 

The time for convection from to the radius r can be found by integration of Eqn. 9. 1 and 
is given by 



dr 1 r^-r'^ 1 



Jrs v{r) IPsC Vs 2psC Vs 

Again, is the expected Alphen velocity of a relativistic plasma. 

Using this equations in Eqn. 9.2, the criterion for a flow to be dominated by convection is 

A± > 3^,r, (9.7) 
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and applying the Steenberg limit of Eqn. 9.5 yields 

PL > 3i8,r,. (9.8) 

Sine is determined through Eqn. 9.6 one obtains by integrating to the radius of the PWN 

(r = %) 



iPscTciRN)' 

Here Tc{Rn) = 1700 y is given by the age of the PWN and the radius of the PWN can be 
estimated from its angular size of 0Ar = 0.1° and its distance d = 5.2 kpc as 

R^ = -^xd=—— x5.2kpc = 9pc. (9.10) 

360° 360° f f 

With these numbers for R]\/ and Tc{Rn) one obtains 

= 1 . = 0.04pc^ll". (9.11) 

2/3sC IVOOy 2x/3(17001y/3.261y)pc 

This is a remarkable result, since it predicts an angular size of 0^ = 11" for the pulsar wind shock 
radius which coincides with feature 5 of the Chandra X-ray image in Fig. 2.5. Therefore this 
theoretical prediction supports an interpretation of the observed ring-like feature in a distance 
of 17" from the pulsar as shock front — similar to the one observed for the Crab Nebula in 
Fig. 2.16. 

Now the right hand side of Eqn. 9.8 provides 

3j8,r, = 0.2pc = 7x lO^^cm. (9.12) 

On the other hand, pi = Eg/eBc is determined by the energy of the electrons (Eg) in the pulsar 

wind where E^, can be estimated from the IC y-ray spectrum measured by H.E.S.S. As discussed 
in the previous section, most of the y radiation can be explained by IC scattering of infrared 
photons with Eph = 1.24 x 10"^ 100 jUm) and CMB photons with Ep^ = 2.35 x lO"'* eV, 

as shown in Fig. 9.1. For these photons the energy of the electrons can be estimated in the 
Thomson limit. According to Eqn. 2.37, one finds 

Ee = meC^{Ey/Ephfl^ (9.13) 

and obtains Eg ~ 2 x (£y/TeV)V2TeV for infrared photons and ~ 20 x {Ey/TdWyl^T€W 
for the CMB photons (cf. Eqn. 2.53). Therefore IC / radiation in the energy range from 0.2 and 
20TeV is mainly produced by electrons of energies between 1 to lOTeV and 10 to 100 TeV (cf. 
Aharonian et al. [2006b]). So in both cases. Eg < 100 TeV and the gyroradius of the electrons 
can be estimated as 

Ee lOOTeV fEyVl^fll^G\ ^ , 

Since the gyroradius of the electrons in the pulsar wind is more than an order of magnitude 
smaller than required by Eqn. 9.8 and Eqn. 9.2 for a diffusion dominated transport, the transport 
cannot be dominated by diffusion. 

This picture of MSH 15—52 is consistent with the general model by de Jager [2006a] ac- 
cording to which PWNs are described as "filled bags" which confine their VHE electron wind 
by the magnetic field which is wrapped in spirals around the pulsar. Due to the perpendicular 
orientation of the magnetic field lines to the radial direction, diffusion of the VHE electrons 
wind is suppressed and a PWN can only release it slowly via convective flows. 
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9.3 Interaction of the Pulsar Wind Nebula with RCW 89 

The sky maps of MSH 15—52 in Sec. 8.6 have provided first insights into the TeV y-ray mor- 
phology of MSH 15—52 which also allow for some further conclusions. 

First, the compact structure extending southeast of the pulsar coincides in size and orien- 
tation with the jet of the pulsar, which has been observed in X-rays. It can therefore can be 
identified as such, which is the first time that a jet has been resolved at TeV energies. 

Second, the region of RCW 89 shows a less intense VHE y-ray emission than the region of 
PSR B 1509— 58, while both regions show a similar intensity in X-rays. A possible explanation 
for this could be that the TeV y radiation is being absorbed in the dense nebula of RCW 89 or 
that the PWN has not expanded far into this region. In the case that the latter is true, one could 
speculate that RCW 89 might constitute a kind of barrier for the PWN that also influences its 
future evolution. The difference between the two regions is illustrated by the multiwavelength 
data of Fig. 9.2 to 9.4. Fig. 9.2 shows a superposition of the 843 MHz Molonglo data (red) 
of the surrounding supernova remnant (Whiteoak and Green [1996]), the ROSAT X-ray con- 
tours (Trussoni et al. [1996]) and the smoothed H.E.S.S. data (green, Aharonian et al. [2005b]). 
Fig. 9.3 shows the Chandra X-ray data (Gaensler et al. [2002]) separated into four different 
energy bands: 6-8 (red), 4-6 (blue), 2-4 (purple) and 0.3-2 keV (green). Fig. 9.4 shows the data 
from the COSMOS H-alpha Survey (Parker et al. [2005]) overlaid with the Chandra X-ray data. 

Third, it is noteworthy that the y-ray measurements are able to indicate a difference between 
the PWN and the optical nebula RCW 89 which is not as obvious from the X-ray measurements 
alone. On the other hand, the resolution of the y-ray data is not alone sufficient enough to 
identify the pulsar jet without the high-resolution images of the X-ray data. In this respect the 
y-ray and X-ray data are complementary to each other. Here multiwavelength data has provided 
new insights, which would have been more difficult to achieve with data from a singe energy 
band only. 

This progress supports hopes that future experiments will make new valuable contributions 
to the understanding of more aspects of MSH 15—52. For example H.E.S.S. II will be able 

to resolve the VHE y-ray morphology at lower energies and with higher resolution than any 
other experiment before. Hopefully these new insights will also increase the knowledge and 
understanding of pulsars and PWNs in general. 
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r ~ ^ 



Ri«ht Ascension (jaOOO) 

Figure 9.2: 843 MHz Molon- 
glo data (red), ROSAT X-ray 
contours and smoothed H.E.S.S. 
data (green). (Figure taken from 
Gaensler and Slane [2006].) 



Figure 9.3: Chandra X-ray 
data at different energy 
bands. (Figure taken from 
NAS A/CXC/MIT/B .Gaensler 
et al. [2006].) 



Figure 9.4: Image from 
the COSMOS H-alpha Survey 
overlaid with Chandra X-ray 
contours. (Figure taken from 
Yatsu et al. [2006].) 
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Chapter 10 



Summary 



This work has reported on the discovery of VHE 7 radiation from the pulsar wind nebula MSH 
15—52 with H.E.S.S. MSH 15—52 is also known as supernova remnant G320.4— 1.2. It hosts 
and is powered by PSR B 1509— 58, one of the most energetic pulsars known, which is embed- 
ded in the supernova remnant. As a prime target, MSH 15—52 was one of the first objects to be 
observed with H.E.S.S., which is currently one of the most sensitive instruments for VHE 7-ray 
astronomy. 

A review of the available data of MSH 15—52 and PSR B1509— 58, since its discovery in 
1961 to the most recent observations with X-ray satellites like Chandra and INTEGRAL shows, 
that MSH 15— 52 is a complex and only partially understood system, which poses many open 
questions to astrophysicists. 

Imaging atmospheric Cherenkov technique provides means for conducting 7-ray astron- 
omy by allowing for the reconstruction of direction and energy of VHE 7-ray photons using 
Cherenkov light which is produced when these photons enter the atmosphere. The detection 
of Cherenkov light is achieved with telescopes with large reflecting dishes. Combining mul- 
tiple telescopes for stereoscopic observations provides further advantages — most importantly 
increased accuracy and sensitivity. 

H.E.S.S. is a stereoscopic system of imaging atmospheric Cherenkov telescopes which has 
been developed to have an unprecedented sensitivity and precision. It is located in Namibia, 
where it has an ideal view of the galactic center. It consists of four telescopes, each with a mirror 
surface of 107 m^, a height of 28 m and a camera with 960 photo-multiplier tubes. H.E.S.S. is 
sensitive in the energy range from 0.2 to 50 TeV and able to detect a 7-ray point source with a 
flux of 1% of the flux from the Crab Nebula in about 25 hours or a source of similar strength as 
the Crab Nebula within 30 seconds. H.E.S.S. has been in operation since June 2002. 

Several standard methods of the imaging atmospheric Cherenkov technique are used in 
H.E.S.S. data analysis and have been introduced here. They include the production of Monte 
Carlo simulations, the reconstruction of 7-ray photons, statistical methods for the calculation 
of significance, position and size of a source and spectroscopy. The methods have been tested 
with data from the Crab Nebula and provide reliable results. Image deconvolution using the 
Richardson-Lucy algorithm has been developed with Monte Carlo simulations and applied here 
to H.E.S.S. data for the first time. In particular for strong sources, interesting results with an 
seemingly higher resolution are obtained. The reconstruction of light curves from pulsars with 
ephemeris from the archive of the Australia Telescope National Facility has been discussed and 
was verified with optical data from the Crab Nebula. 
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CHAPTER 10 Summary 



The methods of data analysis have been applied to the first H.E.S.S. data of MSH 15—52. 
The observations were conducted in 2004 from March 26 to July 20. After run selections, data 
with a live-time of 26.14h was available. MSH 15-52 was detected as HESS J1514-591 with 
a significance up to 32 standard deviation, corresponding to 6 cr/h and a y-ray rate of (2.6 ± 
0.1) min~^ The center of the emission region was found at (J2000) (15hl4'^6.^5 ± 2.Hstat ± 
2.^6syst, -59° 10' l'/2 ±2i;;^t ±20'/y,t) at a distance of 2.3' from PSR B 1509-58. The y-ray 
excess can be fit by a Gaussian distribution with a width and length of 2.3' ± 0.4'^^^^ and 6.5' ± 
0.5stat respectively. The major axis of the emission region has an angle of 43° ± 4° to the RA- 
axis and is aligned with the jet axis of PSR B 1509— 58, which has been observed in X-rays. 

The energy spectrum obeys a power law of the form ^ = ^ixeV (t^v) with ^ixeV = 
(5.8±0.2stat± 1.3syst) X 10-i2cm-2s-iTeV-i, T = 2.32±0.04stat±0.10syst and *(E > ITeV) 
= (4.4 ± 0.2stat ± 1 -Osyst) X 10 ^^cm^s ^ . If the integrated flux above 1 TeV is compared with 
it, MSH 15—52 emits 20% of the flux of the Crab Nebula with a harder photon index. No 
spatial variation of the photon index nor temporal variability of the flux on the timescales of 
days within the observation period was found. 

The Richardson-Lucy algorithm for image deconvolution has been applied to the H.E.S.S. 
y-ray maps and images have been obtained which are complementary to and in agreement with 
those obtained by the standard smoothing technique. They show the VHE y-ray morphology of 
MSH 15—52. A study of the energy dependence of the morphology shows a decreasing y-ray 
emission along the pulsar jet axis with increasing energy. The correlation between the y- and 
X-ray emission was studied based on the X-ray data from the ROSAT and the Chandra satellite. 
While a correlation between y- and X-ray data was found in the region of PSR B 1509— 58, a 
correlation was not found in the northwest region of the optical nebula RCW 89. Between the 
H.E.S.S. and Chandra data, the correlation coefficients of 0.58^q |y and -O.OGto^t 

have been 

obtained for the regions of PSR B 1509 58 and RCW 89, respectively. 

The y-ray excess was tested for pulsed emission from PSR B1509— 58. Since no indication 
for a pulsed emission were found, upper limits for the pulsed y-ray flux were determined for a 
confidence level of 99%. They are 8.3 x 10~^^cm~^s~^ and 0.82x 10~^^cm~^s~^ for the energy 
range above 280 GeV and 1 TeV respectively. 

The results obtained from the data analysis allow various conclusions to be drawn about the 
PWN. For example, it is possible to reproduce the multiwavelength spectrum from MSH 15—52 
with a population of VHE electrons, where the y-ray energies are produced by IC scattering of 
photons from the interstellar radiation field. A hadronic production scenario for the y radiation 
by nuclei of the PWN cannot be excluded, but no evidence for this assumption has been found. 

Moreover, following the argumentation by de Jager [2006c] which is based on the model 
by Kennel and Coroniti [1984], it is possible to show that the transport in the PWN must be 
dominated by convection rather than diffusion. 

Finally, using the correlation between the y-ray and X-ray data, it is possible for the first 
time to detect a jet in TeV y radiation. In addition, the missing correlation in the region of 
the dense optical nebula RCW 89 allows for speculation as to whether or not IC radiation is 
produced in this region and absorbed by the dense nebula of RCW 89, or if instead the PWN 
has not expanded far into this region and RCW 89 constitutes a kind of barrier for the PWN. 
Further experiments such as H.E.S.S. II are expected to provide new data and new information 
about the y-ray morphology of MSH 15-52 as well as about PSR B1509-58 and PWNs in 
general. 
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Appendix A 



Flowchart of the H.E.S.S. Standard 
Analysis Chain 



r SHESSCONFIG 
wdl/hess/analysis/desert 

Effective Areas 

EffectiveAreas.root 



Lookup tables 

energy_info.root 



nCT \*-4 SHESSDSTDATA 
I /dl/hess/dstdata 

T 

hdanalysis 
wobble_chain 

f 




Significance 
map 

Spectrum 
Light curve 



Figure A.l: Flowchart of the H.E.S.S. standard analysis chain showing the dependencies and the 
order of execution of the individual steps. Further details can be found in Klages et al. [2005]. 
(Flowchart taken from Klages et al. [2005].) 
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Appendix B 
Hillas Parameters 



The Hillas parameters were first introduced by Hillas [1985]. They are derived from the mo- 
ments of the image amplitude distribution. If Z, is the amplitude of the pixel / (/ e N) with the 
coordinates jc,- and yt, then the nth order moment (jc^_y^) is defined as 

W) = ^^H^, (B.l) 
with n — p + q and n,p,q^ N. Typical Hillas parameters are shown in Fig. 5.3. 



B.l 0th Order Moment 

The 0th order moment is the image amplitude {I A) given as 

M = £/,. (B.2) 

The image amplitude or size represents the total intensity of the image. 

B.2 1st Order Moments 

The 1st order moments are 

(;c) = ^A, (B.3) 

{y) = (B.4) 

These are the coordinates of the center of gravity. The local distance {LD) is the distance from 
the center of gravity to the camera's origin. It can be expressed using the 1st moments as 

LD = ^(x)2+(3;)2. (B.5) 
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B.3 2nd Order Moments 



APPENDIX B Hillas Parameters 



B.3 2nd Order Moments 



The 2nd order moments are 



which provide the image distributions 



Jyy 



■'xy 



and the auxiliary parameters 



A, 



{xy)-{x){yy 



n = ^ m-^ + 4(xy)2, 

M± = l±k/l. 



(B.6 
(B.7 
(B.8 



(B.9 
(B.IO 
(B.ll 



(B.12 

(B.13 

(B.14 

(B.15 
(B.16 



Then the width (W) and length (L) of the image with respect to the major and minor axes of the 
Hillas ellipse of the image can be written as 



L = ^(0,2 + 0^2 +0/2. 
Finally, the miss, the angles a and and the azwidth are given by 



miss 



= tan 



^J {u+ {x)^ + u- {y)^)/2 - 2{xy)G^/lo, 
-1 f{k + l){y)+2o,y{xy 



a = 



sm 



2a^^{y)-(k-l){x) 
.1 / miss 



-D 



azwidth = {{x)^ + (y)^ —no)/2. 



(B.17 
(B.18 

(B.19 
(B.20 

(B.21 
(B.22 



134 



Appendix C 
Gaussian Fit Function 



The position and size of a y-ray signal is determined from a fit of a two-dimensional Gaussian 
function G to the excess map of the data. G is obtained by a convolution of the two-dimensional 
Gaussian function goj^Gy with the H.E.S.S. point spread function (PSF) to allow a separation 
between the influence of the PSF and true or intrinsic width and length of the source. The 
derivation of G is given here. 

If ^ is a function describing a Gaussian distribution with a zero-mean, then a one-dimen- 
sional representation with standard deviation Ox is 

gaM = ^S^exp (-1-^] . (C.l) 



V27tOx \ 2 C7/ 

The two-dimensional extension is 

ga,,ay {x, y) = go. (x) • gay {y) = ^^^^ e^P f - o I + ^ I I • (^.2) 

Furthermore, with * denoting the convolution operator, the convolution of two Gaussian func- 
tions is 

1 f \ \ 

8oA^)*8c,{x) = :exp ---——J = „ ^^(^) (C.3) 

and consequently 
go„Oy{x,y)*gai,ai{x,y) = [goAx) ■ gay{y)]*[gcji{x) ■ ga,{y)] 

= [<?a, W * gcTi (x)] ■ [goyiy) * -go, iy)] 




1 V y 1 

According to Sec. 5.2.6, the H.E.S.S. PSF (PSF) can be described by the sum of two radial 
symmetric Gaussian functions with normalization constants A and A^ei as 

PSF{x,y) = A[ga, .tTl {X, 

y) + y)] (C.4) 

= Nigai,oi{x,y)+N2ga2,02ix,y) 

Ni f ^ y^W N2 { 1 { x^ / 
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The convolution of gGx,(7y{x,y) with PSF is 



Ga,,ayix,y) = ga,,ayix,y)*PSF{x,y) 

= ga^,ay{x,y) * [Mg<7i,(Ti ix,y) +N2ga2,G2ix,y)] 



Ni 



1 \/ 

N2 



:exp 



2;r A /a2 + a| 



a 2 + a I 



exp 



iV2 



:exp 



:exp - 




+ ■ 



+ ■ 



((72 + (7|) ((72 + (7|) 



+ ■ 



(a2 + C72) (C72 + C72) 



+ ■ 



(a2 + a2) (a2 + a2) 



(C.5) 



According to Eqn. C.5 the relative amplitude A ^g/ before the convolution is given as 



(C.6) 



Consequently, by substituting A^i/A^2 with Eqn. C.6 Gc^^Cy{x,y) can be written as 



Gc,,cjy{x,y) 



exp 



+ ■ 



2 1(ct2 + Ct2) (ct2 + (t2) 



+Arel 



(C.7) 



<^<7r,<7j,(-^,y) is centered at the origin, and the major and minor axes are parallel to the coordinate 
axes. To fit excess distributions of arbitrary orientation, a coordinate transformation by the 

rotation matrix M^o provides the necessary rotation by the angle ft). Additionally, two parameters 
(ao and 5o) are required for fitting an excess with an arbitrary position. The final fit function 
Gax,ay,a(),do,(o,N{^-i ^) equatorial coordinates a and 5 is obtained from substitution of x and y 
as 



a-Oo 
d-do 



cos(ftj)cos(5)(a - oco) - sm{co){d - do) 
sin(ft)) cos(5) (a - oib) +cos(ft))(5 - do) 



(C.8) 
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Appendix D 

Upper Limits According to Feldman and 
Cousins 



The unified approacli by Feldman and Cousins [1998] is a method for determining the confi- 
dence intervals (CI) of a measurement for a given confidence level (CL). A discussion of this 
approach is also found in the Particle Data Booklet (Yao et al. [2006], Chp. 31). The Feldman 
Cousins approach has two important properties: it provides the correct "coverage" and it avoids 
"flip-flopping". Correct coverage means that a measured quantity x lies within the CI with the 
same probability that is given by the corresponding CL, if the measurement is repeated numer- 
ous times. Flip-flopping refers to a problem that arises if different methods for the calculation 
of the CIs of upper limits and measurements are applied, which results in an incorrect coverage. 
The avoidance of flip-flopping is of particular importance for measurements of small signals 
where it is not clear beforehand whether one will find a signal or an upper limit. Flip-flopping 
is avoided in the unified approach by Feldman and Cousins, which provides a unified construc- 
tion of CI regardless of whether a signal or an upper limit is obtained. Here the problem is 
illustrated and the approach is discussed in the case of a "bounded Gaussian". The same dis- 
cussion is found in greater detail in e.g. Feldman and Cousins [1998] and Schwanke and Lohse 
[2004]. 

The classic approach constructs the CI for a CL of a by determining the parameter L such 
that 

n, X \ fnlr^ Pi-x\u)dx = a, in the case of a signal , 
I J_^P{x\iJ,)dx = a, in the case of an upper limit. 

However, the transition from one case to the other implies flip-flopping and does not guarantee 
correct coverage. 

Therefore, the approach by Cousins and Feldman introduces an ordering parameter R, which 
specifies in what order the CI is constructed. R is defined through the Likelihood ratio of the 
experiment's probability density function (PDF) P{x\il) and P{x\jl) as 

« = ^^, (D.2) 

where x and ji are the measured and true values and fi denotes the value of jU which maximizes 
P(x|ju) and which is physically allowed. Acceptance intervals [xi,X2] are then constructed for 
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= 1 




6 8 10 



Figure D.l: Ordering parameter R{x) for a 
bounded Gaussian PDF (Eqn. D.6). The red, 
green and blue curves correspond to n =0.1, 
0.5 and 5.0, respectively. (Figure taken from 
Schwanke and Lohse [2004].) 




4 6 8 
Measured Mean x 



Figure D.2: Confidence belts for a bounded 
Gaussian (Eqn. D.6) with a CL of 90% con- 
stituted by a set of acceptance intervals. The 
blue confidence belt includes systematic errors; 
the red does not. The corresponding CI [Aii,At2] 
is determined through the intersections with 
the vertical line representing a measured value 
X. (Figure taken from Schwanke and Lohse 
[2004].) 



all possible measurements of x to fulfill 

P{x e [xi,X2]\n) = / P{x\n)dx = a, 
and including the additional conditions 

R{xi) =R{X2). 



(D.3) 



(D.4) 



The union of all acceptance intervals of the possible values of x will provide a so-called confi- 
dence belt as shown in Fig. D.2 from which the CI [/li,jU2] is determined. 

In the case of a Gaussian PDF where the physical values of jU are restricted to positive values 



P{x\il) 



1 



exp 



(D.5) 



For a measurement x, P{x\fl) is maximized ii fi = x. This is true for positive values of x, 
including zero. If x is negative, /i = is the physically allowed value which maximizes P{x\fl), 
i.e. 



R{x) 



exp(-^(x-/i)^), .x:>0 
exp(-i(x-/i)^ -jc^), x<0. 



(D.6) 



The corresponding values of R for =0.1, 0.5 and 5.0 are shown by the red, green and blue 
curves in Fig. D.l. For values of /i that are small in comparison to the standard deviation 
a = 1, the asymmetry is apparent, which is reflected in the acceptance interval. The resulting 
confidence belt for a CL of 90% is shown in Fig. D.2 by the red curves. The CI [/ii,/i2] is 
determined by the intersection of the vertical line at the measured value x with the confidence 
belt, as shown in Fig. D.2. If the vertical line at x does not intersect the lower boundary at 
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Hi > 0, one does not obtain a lower boundary but only an upper limit. The Feldman Cousins 
approach therefore guarantees a smooth transition between a measurement and an upper limit 
and it avoids any discontinuities and flip-flopping. 

It is also possible to incorporate systematic errors into this approach. Therefore it is assumed 
that the systematic error (cTy) has a Gaussian PDF 

where s — (1 ±e) represents a factor by which an actual measurement could appear multiplied 
and where e is of the order of a^. With this ansatz one can formulate a new PDF 

F{s\Gs) = J P{x\sfl)L{s\Gs)ds (D.8) 

which replaces the PDF in Eqn. D.5. For the correct treatment one also has to replace R accord- 
ingly, i.e. 

With these modifications, one proceeds as described above. The CI with the same parameters 
but with an additional systematic error of 20% is shown by the blue curves in Fig. D.2. The new 
confidence interval is larger. 

In the H.E.S.S. standard analysis, the calculation of the CI for a bounded Gaussian according 
to the unified approach by Feldman and Cousins is realized by the class TBoundedGaussian. 
This class is initialized with a CL a and a systematic error Og. The CI is obtained through the 
method GetConfidenceInterval(x, Ox, , &112), which returns the boundaries of the CI jUi and 
II2 for the measured value x and the error Ox- 
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Appendix E 



Lists of Observation Runs 



Table E.l: Run list of H.E.S.S. wobble observation at (J2000) (15^14^27% -59°16'18") of the Crab 
Nebula. The observation period lasted from January 25, 2004 to March 4, 2005. 





Run no. 


Wobble onset 


<Amiuae>[ J 


<Kate>[rlzj 


Events 


Duration[s] 


1 


18526 


^ ,^ A CO 

Dec —0.5 


C A O 

59.3 


o o 1 
231 


A511b5 


1 /COO 

1682 


2 


18870 


Uec +0.5 


4o.Z 


294 


CC/C 1 A/C 

556196 


1 /COO 

1683 


3 


1 oon A 

loo74 


"r\__ 1 (\ CO 

Dec +0.5 


60.1 


OAO 

208 


OA'l O AO 

393848 


1 /CO 1 

1681 


4 


1 oonc 


Dec —0.5 


CA n 


176 


341233 


1 /COO 

1683 


5 


1 o o on 

18889 


T-v 1 r\ CO 

Dec +0.5 


CO A 

58.0 


O 1 1 

213 


406677 


1 /CO'? 

1683 





18890 


T~\-,^ A CO 

Dec —0.5 


/CO C 

62.5 


IOC 

185 


O C 1 CI o 

351518 


1 /CO ^5 

1683 


1 


18927 


T-x__. 1 A CO 

Dec +0.5 


61.1 


1 AC 

195 


372168 


1 /CO /i 

1684 


o 
o 


23037 


T^„„ A CO 

Dec —0.5 


CO £ 

5o.6 


1 1 '3 

113 


1 AO^A'3 

198793 


1 can 
1567 


n 

y 


23062 


T> A A CO 

KA —0.5 


61.9 


83.3 


1 A nc^Ci 

147520 


1 con 

1687 


1 A 

ID 


23063 


DA 1 A CO 

RA +0.5 


cn n 


1 A 1 

101 


1 OA'? A A 

180349 


1 /CO/C 

1686 


1 1 


O o o A/1 

23304 


Ty \ A CO 

RA —0.5 


CO A 

52.4 


137 


O /I O AC /I 

243054 


1 /C 1 o 

1613 






uec -rU.j 


'to. J 


1 89 


a/1 1 CI 


1 ^.86 
iUoO 


13 


23310 


Dec -0.5° 


50.2 


171 


323219 


1686 


14 


23526 


Dec +0.5° 


46.3 


213 


425821 


1684 


15 


23544 


RA +0.5° 


48.8 


210 


435015 


1686 


16 


23545 


RA -0.5° 


46.3 


220 


460561 


1686 


17 


23546 


Dec +0.5° 


46.0 


221 


264014 


948 


18 


23547 


Dec -0.5° 


44.8 


227 


270707 


942 


19 


23577 


RA -0.5° 


45.3 


214 


451231 


1686 


20 


23579 


RA -0.5° 


47.3 


208 


246834 


972 


21 


23593 


Dec -0.5° 


46.2 


204 


402139 


1687 


22 


23642 


Dec -0.5° 


51.1 


193 


379326 


1686 


23 


23662 


Dec +0.5° 


46.0 


218 


3AT2QA 


1679 


24 


23739 


RA +0.5° 


46.2 


189 


373756 


1686 


25 


23740 


Dec -0.5° 


44.9 


192 


381556 


1686 


26 


23741 


Dec +0.5° 


46.2 


186 


372221 


1686 


27 


23753 


RA +0.5° 


52.1 


169 


316510 


1686 


28 


23756 


Dec +0.5° 


46.8 


199 


408124 


1687 


29 


24138 


Dec +0.5° 


48.1 


192 


187203 


907 


30 


24412 


RA +0.5° 


51.8 


166 


308661 


1687 


30 






51.5 


190 


10326374 


47385 
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Table E.2: Run list of all H.E.S.S. observation runs of MSH 15—52 from the observation period of 
March 26 to July 20, 2004 taken at the wobble position (J2000) (15''14"'27^ -59°16'18") +0.5° in 
Dec. Runs in parentheses did not pass the run selection criteria and were not used in the analysis. 





Run no. 


<Altitude> [°] 


<Rate> [Hz] 


Events 


Duration [s] 


- 


(20136) 


(51.9) 


(287) 


(520570) 


(1682) 


1 


20282 


53.4 


259 


480314 


1687 


2 


20301 


53.7 


269 


502145 


1683 


3 


20303 


54.3 


272 


521903 


1682 


4 


20323 


54.4 


295 


586611 


1682 


5 


20325 


53.5 


294 


580954 


1682 


6 


20343 


53.2 


286 


540909 


1682 


7 


20345 


54.4 


289 


550983 


1683 


8 


20366 


54.4 


261 


487784 


1683 


9 


20368 


53.3 


258 


485038 


1683 


10 


21085 


53.1 


249 


547372 


1683 


- 


(21181) 


(54.0) 


(246) 


(457715) 


(1684) 


11 


21207 


53.3 


233 


475160 


1683 


12 


21210 


54.4 


237 


490840 


1683 


13 


21228 


45.2 


199 


375716 


1682 


- 


(21230) 


- 


- 


(0) 


(1681) 


- 


(21232) 


- 


- 


(0) 


(1681) 


14 


21261 


51.2 


227 


425355 


1683 


15 


21263 


54.1 


239 


446111 


1683 


16 


21266 


53.8 


238 


443065 


1683 


17 


21290 


54.4 


227 


419622 


1683 


18 


21368 


54.2 


223 


412307 


1683 


19 


21580 


54.4 


222 


436539 


1686 


20 


21598 


53.9 


227 


476823 


1683 


21 


21618 


54.2 


207 


397984 


1683 




(21641) 


(54.3) 


(96) 


(167231) 


(1684) 


22 


21643 


52.0 


200 


390105 


1686 




(21688) 


(54.5) 


(215) 


(391629) 


(1681) 




(21690) 




(213) 


(54357) 


(1685) 




(21711) 






(0) 


(424) 


23 


21739 


49.8 


191 


350838 


1685 


23 




53.3 


245 


10824478 


40407 
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Table E.3: Run list of all H.E.S.S. observation runs of MSH 15—52 from the observation period of 
March 26 to July 20, 2004 taken at the wobble position (J2000) (1544"^27% -59°16'18") -0.5° in 
Dec. Runs in parentheses did not pass the run selection criteria and were not used in the analysis. 
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Table E.4: Run list of all H.E.S.S. observation runs of MSH 15—52 from the observation period of 
March 26 to July 20, 2004 taken at the wobble position (J2000) (15*'14'"27% -59°16'18") +0.5° in 
RA. 





Run no. 


<Altitude> [°] 


<Rate> [Hz] 


Events 


Duration [s] 


1 


20391 
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283 


542630 


1682 
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1683 



7 - 53.3 277 3659693 11779 



Table E.5: Run list of all H.E.S.S. observation runs of MSH 15—52 from the observation period of 
March 26 to July 20, 2004 taken at the wobble position (J2000) (15'>14'"27% -59°16'18") -0.5° in 
RA. Runs in parentheses did not pass the run selection criteria and were not used in the analysis. 
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Appendix F 

The Richardson-Lucy Algorithm 



Listing F.l: The "max_likelihood.pro" procedure as provided by the "IDL Astronomy Library" 
(Landsman [2004]) for the Richardson-Lucy deconvolution. 

+ 

NAME: 

MAX_L IKEL I HOOD 
PURPOSE: 

Maximum likelihood deconvolution of an image or a spectrum. 
EXPLANATION: 

Deconvolution of an observed image (or spectrum) given the 
instrument point spread response function (spatially invariant psf) . 
Performs iteration based on the Maximum Likelihood solution for 
the restoration of a blurred image (or spectrum) with additive noise. 
Maximum Likelihood formulation can assume Poisson noise statistics 
or Gaussian additive noise, yielding two types of iteration. 

CALLING SEQUENCE: 

for i=l, Niter do Max_Likelihood, data, psf, deconv, FT_PSF=psf_ft 

INPUTS PARAMETERS: 

data = observed image or spectrum, should be mostly positive, 

with mean sky (background) near zero, 
psf = Point Spread Function of the observing instrument, 

(response to a point source, must sum to unity) . 

INPUT/OUTPUT PARAMETERS : 

deconv = as input: the result of previous call to Max_Likelihood, 
(initial guess on first call, default = average of data), 
as output: result of one more iteration by Max_Likelihood. 

Re_conv = (optional) the current deconv image reconvolved with PSF 
for use in next iteration and to check convergence . 



OPTIONAL INPUT KEYWORDS: 

/GAUSSIAN causes max-likelihood iteration for Gaussian additive noise 
to be used, otherwise the default is Poisson statistics . 
FT_PSF = passes (out /in) the Fourier transform of the PSF, 

so that it can be reused for the next time procedure is called, 
/NO_FT overrides the use of FFT, using the IDL function convol () instead. 
POSITIVITY_EPS ^ value of epsilon passed to function posltivity, 

default = -1 which means no action (identity) . 
UNDERFLOW_ZERO = cutoff to consider as zero, if numbers less than this. 

EXTERNAL CALLS: 

function convolve ( image, psf ) for convolutions using FFT or otherwise, 
function posltivity ( image, EPS= ) to make image positive. 

METHOD: 

Maximum Likelihood solution is a fixed point of an iterative eq. 
(derived by setting partial derivatives of Log (Likelihood) to zero). 
Poisson noise case was derived by Richardson (1972) & Lucy (1974) . 
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; Gaussian noise case is similar with subtraction instead of division. 

; HISTORY: 

; written: Frank Varosi at NASA/GSFC, 1992. 

; F.V. 1993, added optional arg. Re_conv (to avoid doing it twice). 

; Converted to IDL V5 . W. Landsman September 1997 

r 

pro Max_Likelihood, data, psf, deconv, Re_conv, FT_PSF=psf_ft, NO_FT=noft, $ 

GAUSSIAN=gaussian, $ 
POSITIVITY_EPS=epsilon, $ 
UNDERFLOW_ZERO=under 

if N_elements ( deconv ) NE N_elements( data ) then begin 
deconv = data 

deconv [*] = total ( data ) /N_elements ( data ) 
Re_conv = 

endlf 



if N_elements( under ) NE 1 then under = l.e-22 
if N_elements ( epsilon ) NE 1 then epsilon = -1 

if N_elements( Re_conv ) NE N_elements ( deconv ) then $ 

Re_conv = convolve ( positivity( deconv, EPS=epsilon ), psf, $ 

FT_PSF=psf_ft, NO_FT=noft ) 



if keyword_set ( gaussian ) then begin 



deconv = deconv + convolve ( data - Re_conv, psf, /CORREL, $ 

FT_PSF=psf_ft, NO_FT=noft ) 

endif else begin 

wp = where ( Re_conv GT under, npos) 
wz = where ( Re_conv LE under, nneg) 

if (npos GT 0) then Re_conv[wp] = ( data [wp] /Re_conv [wp] ) > 
if (nneg GT 0) then Re_conv[wz] = 1 

deconv = deconv * convolve ( Re_conv, psf, FT_PSF=psf_ft, $ 

/CORREL, NO_FT=noft ) 

endelse 



if N_params ( ) GE 4 then $ 

Re_conv = convolve ( positivity{ deconv, EPS=epsilon ), psf, $ 

FT_PSF = psf_ft, NO_FT = noft ) 

end 
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Appendix G 

PSF and Source Extension at Different 
Energy Bands 



Tbl. G.l shows the individual best fit parameters (cf. Eqn. 5.11) and the containment radius 
(''68%) of the H.E.S.S. PSF for different energy bands. They have been determined from point 
source Monte Carlo data with zenith angle of 40°, wobble offset of 0.5° and a restriction to 
the corresponding energy range. The variation of the PSF is explained by the reconstruction 
accuracy for y-ray showers which increases with energy. The parameter Os denotes the standard 
deviation of the Gaussian which has been applied to smooth two maps to a similar PSF with 
same m%. 

Fig. G.l shows slices along the major and minor axis of the excess maps and the fit function 
G (cf. App. C) of Fig. 8.17 (left column). In addition to G (black) its two components of the 
PSF (red) and the intrinsic width (green) are shown. 



Table G.l: PSFs at different energy bands shown in Fig. 8.17. The parameters A^/, (Ji and (72 are 
determined by Eqn. 5.11. r^^% denotes the 68% containment radius and CT, the standard deviation of 

the smoothing Gaussian appUed before the parameters were determined. The values of the last line 
in brackets show the PSF standard parameters of the unrestricted energy range for comparison. 



E [TeV] Gs Arel gl[°] 02[°] . 

0.2TeV<£:<0.5TeV 0.35 0.0574 0.108 0.127 

0.5 TeV< ^ < 1 TeV 0.17 0.0479 0.110 0.110 

lTeV<£<2TeV 0.096 0.0409 0.116 0.102 

2TeV<£<5TeV 0.057 0.0361 0.123 0.095 

5TeV<£:< lOOTeV 0.052 0.0275 0.0929 0.067 

0.2TeV<£:<5TeV 0.181 0.0461 0.113 0.118 

0.2TeV<£<5TeV 0.023 0.219 0.0518 0.116 0.125 

5TeV<£ < lOOTeV 0.064 0.298 0.0673 0.106 0.125 

(0.2TeV<£ < lOOTeV 0.170 0.0477 0.117 0.120) 
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Figure G.l: Fit functions along the 
major axis (upper plot) and minor axis 
(lower plot) of the y-ray excess maps of 
Fig. 8.17 (left column). Each fit function 
as well as their components of the intrin- 
sic width and the PSF are shown. 
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Appendix H 

Chandra Data Analysis 



The Chandra X-ray data discussed in Sec. 8.6.3 was obtained from the Chandra Data Archive 
NASA/CXC [2006]. The data was taken by the four ACIS CCD detectors on 2000 August 14 
(observation ID 754) during a single exposure of 20 ks. After standard processing at the Chandra 
Science Center, the further analysis was carried out with the "Chandra Interactive Analysis of 
Observations" (CIAO 3.3, Fruscione et al. [2006]). The list of bad pixels was provided with the 
data and taken into account in the analysis. After dead-time corrections a live-time of 19039 s 
with a total of 259770 events remained. A further restriction of the energy range to 0.3-lOkeV, 
where the detector response is reliable, reduced the data set to a total of 222914 events. The 
corresponding count map is shown in Fig. H.l. The exposure corrected count map of Fig. H.3 
was obtained by a bin-wise devision of the count map by the exposure map which is shown in 
Fig. H.2. The exposure map for this observation period was created following the guidelines of 
the Chandra analysis software. The exposure map clearly shows the four individual ACIS CCD 
detectors separated by gaps of reduced exposure. While these gaps are clearly seen in the count 
map, they have been compensated in the exposure corrected map. 




Figure H.1: Chandra X-ray 
count map of the ACIS CCD 
detectors of MSH 15-52 af- 
ter calibration and dead-time 
correction. 



Figure H.2: Chandra expo- 
sure map for the ACIS CCD 
corresponding to the obser- 
vation of Fig. H.l. 
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Figure H.3: The corre- 
sponding Chandra exposure 
corrected X-ray map of 
Fig. H. 1 . The peak intensity 
is found at the position of 
PSR B 1509-58 which is 
indicated by the black circle. 
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CCD 


charge-coupled device 




confidence level 




cosmic microwave background 


i_ 1 


i^nerenKOV icicscopc 




declination 


Hnf 
UOi 


degrees oi irceuum 


FDV 
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TT C C 


High Energy Stereoscopic System 


TA 
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TAPT 
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T THr» 


iigni emiiiing uiocie 


MHD 


magneto-hydrodynamic 


MJD 


Modified Julian Date 


MRSL 


mean reduced scaled length 


MRSW 


mean reduced scaled width 


PMT 


photo multiplier 


PWN 


pulsar wind nebula 


PDF 


probability density function 


PSF 


point spread function 


RA 


right ascension 


RL 


Richardson-Lucy 


RMS 


root mean square 


SNR 


supernova remnant 


VHE 


very high energy (10GeV< E <100TeV) 
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